Testing of Transmitters for Communication Links by Software Simulation of Reference Channel and/or Reference Receiver

ABSTRACT

A transmitter for a communications link is tested by using a (software) simulation of a reference channel and/or a reference receiver to test the transmitter. In one embodiment for optical fiber communications links, a data test pattern is applied to the transmitter under test and the resulting optical output is captured, for example by a sampling oscilloscope. The captured waveform is subsequently processed by the software simulation, in order to simulate propagation of the optical signal through the reference channel and/or reference receiver. A performance metric for the transmitter is calculated based on the processed waveform.

CROSS-REFERENCE TO RELATED APPLICATION(S)

This application is a continuation of U.S. application Ser. No.13/367,234 filed Feb. 6, 2012 entitled “Testing of Transmitters forCommunication Links by Software Simulation of Reference Channel and/orReference Receiver,” which is a continuation of U.S. application Ser.No. 12/430,892 filed Apr. 27, 2009, entitled “Testing of Transmittersfor Communication Links by Software Simulation of Reference ChannelAnd/or Reference Receiver,” which claims priority under 35 U.S.C.§119(e) to U.S. Provisional Patent Application No. 61/048,063,“Alignment of Square Waveform in T/WDP Code for a More Accurate Measureof xMA,” filed Apr. 25, 2008. U.S. application Ser. No. 12/430,892 isalso a continuation-in-part of U.S. patent application Ser. No.11/316,115, “Testing of Transmitters for Communication Links by SoftwareSimulation of Reference Channel and/or Reference Receiver,” filed onDec. 21, 2005 which claims priority under 35 U.S.C. §119(e) to each ofthe following U.S. Provisional Patent Applications Ser. No. 60/638,788,“Proposed TP-2 Test Methodology,” filed Dec. 22, 2004; Ser. No.60/641,834, “Proposed TP-2 Test Methodology,” filed Jan. 5, 2005; Ser.No. 60/643,234, “Proposed TP-2 Test Methodology with MATLAB Script,”filed Jan. 11, 2005; Ser. No. 60/677,911, “New Approach To Measure TxSignal Strength And Penalty,” filed May 4, 2005; Ser. No. 60/690,190,“Improvements To Transmit Waveform And Dispersion Penalty Algorithm,”filed Jun. 13, 2005; Ser. No. 60/698,435, “Method To Accurately MeasureBias And Optical Modulation Amplitude For a Distorted Optical Waveform,”filed Jul. 11, 2005; Ser. No. 60/701,637, “Method To ExtractSquare-Wave-Based Bias And Optical Modulation Amplitude For A DistortedOptical Waveform,” filed Jul. 21, 2005; Ser. No. 60/707,282, “Method ToExtract Square-Wave-Based Bias And Optical Modulation Amplitude For ADistorted Optical Waveform,” filed Aug. 10, 2005; Ser. No. 60/710,512,“Computation Of Optimal Offset For A Transmitter Waveform And DispersionPenalty,” filed Aug. 22, 2005; and Ser. No. 60/713,867, “TransmitterWaveform And Dispersion Penalty Test Using A Short Equalizer, OMSDNormalization, And Optimal Offset Computation,” filed Sep. 1, 2005. Thesubject matter of all of the foregoing are incorporated herein byreference in their entirety.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates generally to the use of software simulation(which is intended to include firmware simulation) in testing oftransmitters used in communications links, for example the testing oftransmitters used in optical fiber communications link by softwaresimulation of a reference channel (i.e., a reference optical fiber link)and/or a reference receiver.

2. Description of the Related Art

Optical fiber is widely used as a communications medium in high speeddigital networks, including local area networks (LANs), storage areanetworks (SANs), and wide area networks (WANs). There has been a trendin optical networking towards ever-increasing data rates. While 100 Mbpswas once considered extremely fast for enterprise networking, attentionhas recently shifted to 10 Gbps, 100 times faster. As used in thisapplication, 10 Gigabit (abbreviated as 10 G or 10 Gbps) systems areunderstood to include optical fiber communication systems that have datarates or line rates (i.e., bit rates including overhead) ofapproximately 10 Gigabits per second.

Regardless of the specific data rate, application or architecture,communications links (including optical fiber communications links)invariably include a transmitter, a channel and a receiver. In anoptical fiber communications link, the transmitter typically convertsthe digital data to be sent to an optical form suitable for transmissionover the channel (i.e., the optical fiber). The optical signal istransported from the transmitter to the receiver over the channel,possibly suffering channel impairments along the way, and the receiverthen recovers the digital data from the received optical signal.

For example, a typical 10 G optical fiber communications link 100 isshown in FIG. 1. The link 100 includes a transmitter 105 coupled throughoptical fiber 110 (the channel) to a receiver 120. A typical transmitter105 may include a serializer, or parallel/serial converter (P/S) 106 forreceiving 10 G data from a data source on a plurality of parallel linesand providing serial data to a 10 G laser driver 107. The driver 107then drives a 10 G laser source 108, which launches the optical waveformcarrying the digital data on optical fiber 110.

On the receive side, a typical receiver 120 includes a 10 Gphotodetector 111 for receiving and detecting data from the opticalfiber 110. The detected data is typically processed through a 10 Gtransimpedance amplifier 112, a 10 G limiting amplifier 113, and a 10 Gclock and data recovery unit 114. The recovered data may then be placedon a parallel data interface through a serial/parallel converter (S/P)115.

Standards play an important role in networking and communications. Sincecomponents in the network may come from different vendors, standardsensure that different components will interoperate with each other andthat overall system performance metrics can be achieved even whencomponents are sourced from different vendors. For example, standardsfor transmitters can be used to ensure that, when compliant transmittersare combined with a compliant channel and compliant receivers, theoverall link will meet certain minimum performance levels. As a result,manufacturers of transmitters typically will want to test theirtransmitters for compliance with the relevant standards as part of theirproduction process.

However, current testing approaches, which were developed for earliergeneration optical fiber communications links, may not be suitable forthe more advanced systems that are currently being developed andfielded. In one testing approach, transmitters are tested based on thequality of their immediate output. For example, the signal to noiseratio (SNR), bit error rate, eye diagram or other performance metric maybe measured at the output of the transmitter. This approach ignoresinteractions between the transmitter and other components in the system(e.g., the optical fiber and the receiver) and, as a result, may be toosimplistic to give an accurate assessment of the overall systemperformance achievable by the transmitter. For example, if the opticalfiber communications link is designed so that the receiver compensatesfor distortions or other impairments introduced by the channel or thetransmitter, a direct measurement of the performance at the transmitteroutput may not accurately reflect this subsequent compensation. In fact,in certain systems, it may be difficult to obtain an accurate assessmentof the performance of the transmitter based only on its immediateoutput. Similarly, standards for transmitters used in more advancedsystems may be based in part on their interaction with other components.For example, if the standard assumes a certain type of distortionmitigation in the receiver, the transmitter may be allowed a certainlevel of distortion that the receiver can compensate, but other types ofdistortion that the receiver cannot compensate should be limited. Atransmitter test that cannot distinguish between the two types ofdistortion (an eye diagram, for example) will be an inadequatecompliance test for this type of standard. Hence, a test of thetransmitter alone may not be sufficient to determine whether thetransmitter complies with the standard.

Another testing approach uses hardware reference components to emulatethe overall optical fiber communications link. For example, there are anumber of standards that relate to 10 G fiber networks. 10 G Ethernetover optical fiber is specified in IEEE Std 802.3ae-2002 Media AccessControl (MAC) Parameters, Physical Layers, and Management Parameters for10 Gb/s Operation (abbreviated herein as IEEE 802.3ae). The IEEE 802.3aestandard includes a Transmitter and Dispersion Penalty (TDP)measurement. The TDP measurement was developed as a transmittercompliance test to ensure that a compliant transmitter will communicatewith a compliant receiver over a compliant channel with acceptableperformance. The TDP measurement involves signal quality characteristicssuch as rise and fall times, eye pattern opening, jitter, and otherdistortion measurements. The intent is to insure that the signalreceived from a compliant transmitter can be detected to deliver aspecified bit error rate.

The TDP measurement technique described in the 802.3ae standard involvesa hardware test setup, including a hardware reference transmitter, thetransmitter under test, a reference channel (i.e., actual length offiber), and a hardware reference receiver. The test proceeds as follows.A reference link is established by connecting the reference transmitterto the reference receiver via the reference channel. The performance ofthe reference link is measured. A test link is established by connectingthe transmitter under test to the reference receiver via the referencechannel. The performance of this link is also measured. The performanceof the reference link is compared with the performance of the test linkto arrive at a value of TDP. The reference transmitter and referencereceiver are high-quality instrument-grade devices. The test procedureinvolves several manual calibration and measurement steps. There areseveral disadvantages to this approach.

One disadvantage is the need for a hardware reference transmitter, ahardware reference channel, and a hardware reference receiver. Thereference transmitter, channel, and receiver have tightly definedcharacteristics and are therefore very expensive. For example, thereference transmitter has tight specifications on rise and fall times,optical eye horizontal and vertical closure, jitter, and relativeintensity noise (RIN). Furthermore, it is possible that differentreference transmitters, channels, or reference receivers might yielddifferent results due to parameter tolerances, environmental conditions(temperature, aging), and other differences. Therefore, the recommendedmeasurement procedure involves a number of complex and sensitivecalibration steps. The accuracy of the TDP measurement is sensitive tothese time consuming and difficult calibration steps. Since thesecalibration steps all have tolerances, there may be significantvariation in the results.

Another disadvantage is that the hardware approach is not easilyextendible to other standards. Another standards committee is IEEE802.3aq, which is developing a standard (10 GBASE-LRM) for 10 G Ethernetover multi-mode fiber over distances of up to 220 meters usingelectronic dispersion compensation (EDC). This standard is in a draftstate, currently documented in IEEE Draft P802.3aq/D3.0, Draft amendmentto: IEEE Standard for Information technology—Telecommunications andinformation exchange between systems—Local and metropolitan areanetworks—Specific requirements, Part 3: Carrier Sense Multiple Accesswith Collision Detection (CSMA/CD) Access Method and Physical LayerSpecifications, Amendment: Physical Layer and Management Parameters for10 Gb/s Operation, Type 10 GBASE-LRM, referred to herein as IEEE802.3aq/D3.0. The use of EDC in the receiver allows communication overlonger distances and/or use of lower cost components. Some of the addedwaveform distortions can be corrected in the EDC receiver. However, thisleads to some problems in defining and measuring a metric fortransmitter qualification.

For example, TDP is not an appropriate transmitter compliance test for10 GBASE-LRM (EDC) systems because the reference receiver does notinclude equalization and does not accurately compensate for thosedistortions that a receiver using EDC can clean up. In general, thereare certain signal impairments that can be corrected by an EDC receiver(such as most linear distortions), and there are other impairments thatcannot normally be equalized (such as most nonlinear distortions). Asuitable transmitter quality measure for an EDC system should measurethe signal quality (or resulting bit error rate) achievable by areceiver with equalization. A reference receiver for such a systemprobably would include a decision feedback equalizer (DFE) as an EDCtechnique. However, a reference grade equalizer in the hardware receiverwould be very difficult to make. It would have its own limitations, suchas finite length and finite precision and tolerances, and theselimitations could conceal the properties of the signal being measured.The resulting measurements may not accurately indicate correctableversus non-correctable impairments.

Another disadvantage with the use of TDP when applied to a transmitterintended for use on a link with an EDC receiver is the dependence of TDPon a measurement of Optical Modulation Amplitude (OMA). OMA is thedifference in optical power for the nominal “1” and “0” levels of theoptical signal. The measurement technique defined by 802.3ae is tocapture samples of a test waveform using a sampling oscilloscope. Thetest waveform is the response of the transmitter to a square wave, whichis a repetitive pattern of several consecutive 1's followed by an equalnumber of consecutive 0's. The mean optical power level of an optical 1is measured over the center 20% of the time interval where the signal ishigh, and similarly for 0's when the signal is low. This is appropriatewhen the waveform is well-behaved in the sense that the levels duringthe measurement windows are constant and well-defined. However, a systemthat uses EDC in the receiver enables the use of lower cost transmittersthat may result in distorted transmitted waveforms, for example asdiscussed further in U.S. patent application Ser. No. 11/033,457, “Useof Low-Speed Components in High-Speed Optical Fiber Transceivers,” filedJan. 10, 2005. The response of the transmitter to a square wave inputmay differ dramatically from a well-behaved waveform due tobandlimiting, dispersion, and other distortions (such as ringing). Thismay make it difficult for an operator to accurately determine the truehigh and low levels. Since the value of TDP is sensitive to the measuredOMA, this inaccuracy will directly affect the TDP measurement.

Thus, there is a need for improved transmitter testing techniques forcommunications links, including for example 10 G optical fibercommunications links where the receiver includes EDC.

SUMMARY OF THE INVENTION

The present invention overcomes the limitations of the prior art byusing a (software—including firmware) simulation of a reference channeland/or a reference receiver to test transmitters. In one embodiment foroptical fiber communications links, a data test pattern (also referredto as a data sequence) is applied to the transmitter under test and theresulting optical output is captured, for example by a samplingoscilloscope. The captured waveform is subsequently processed by thesoftware simulation, in order to simulate propagation of the opticalsignal through the reference channel and/or reference receiver. Aperformance metric for the transmitter is calculated based on theprocessed waveform.

In one embodiment, the performance metric is the transmitter waveformand dispersion penalty (TWDP). TWDP is defined as the difference (in dB)between a reference signal to noise ratio (SNR), and the equivalent SNRat the slicer input of a reference decision feedback equalizer (DFE)receiver for the measured waveform from the transmitter under test afterpropagation through an optical fiber reference channel.

In one approach, a data test pattern is applied to the transmitter undertest and the transmitter output is sampled using standard testequipment. The captured waveform is processed through a softwaresimulation of the reference optical fiber (i.e., reference channel) andthe reference receiver (which includes an equalizer). The use of astandard software simulation of the reference channel and the referencereceiver avoids the difficulties inherent with a hardware receiver andallows a more accurate measurement of unequalizable signal impairments.The need for a hardware reference transmitter, fiber, or referencereceiver is eliminated, thus saving significant costs. In this way, TWDPcan be measured for optical fiber communication systems using EDC.

This approach has advantages since using a hardware reference receiverhas significant problems. For example, the IEEE 802.3aq committeeadopted the definition of one reference receiver based on PIE-D (Penaltyfor Ideal Equalizer—DFE), which corresponds to an infinite length DFEand can be computed analytically. However, an EDC receiver that can bereasonably implemented in hardware will have a finite length equalizerwith limited arithmetic precision. Therefore, the results from ahardware reference receiver will not fully reflect the theoretical TWDPresults and may be strongly influenced by the actual reference receiverimplementation.

In one aspect of the invention, the signal output from the transmitterunder test is captured by a sampling oscilloscope, and the resultingsampled signal is processed by a software simulation (including firmwaresimulations) of a reference channel and reference receiver. This enablesthe TWDP measurement to be performed using standard test equipment,rather than a very expensive instrument grade receiver. In oneimplementation, the reference equalizing receiver is implemented using aMATLAB simulation or other equivalent software simulation. Some or allof the software simulator may be integrated as part of the samplingoscilloscope (or other waveform capture device), provided as an add-onmodule to the sampling oscilloscope and/or implemented as a separatemodule.

The software simulator can also perform measurement of opticalmodulation amplitude (OMA) and waveform baseline. The computation ofthese quantities based on the captured waveform (including earlier orlater processed versions of the captured waveform) eliminates theinaccuracies inherent in separate or manual measurements. For example,these quantities can be computed from the raw digitized waveformcaptured by the waveform capture device, from the waveform after somepre-processing (e.g., after re-sampling, truncation and/or dataalignment) or from waveforms after simulation of propagation. Techniquesare disclosed that yield more accurate measurements for highly distortedsignals. For example, temporal misalignment between a synthesizedwaveform, such as a square wave, and measurement windows capturingvalues for OMA calculation can be corrected, allowing time alignment ofthe synthesized waveform and the measurement windows. This isadvantageous because TWDP depends strongly on OMA, and the accuratemeasurement of OMA may be very difficult for distorted transmitwaveforms that, in spite of this distortion, are acceptable when EDC isused in the receiver. Temporal alignment between the synthesizedwaveform and the measurement windows may also allow more accuratecalculation of waveform and dispersion penalty (WDP) of optical orelectrical signals at different locations within a system. In anothervariation, the software automatically corrects for the optimal DC offset(bias) value. Since the TWDP value is sensitive to offset of thewaveform, automatic offset optimization improves reproducibility of themeasured TWDP. This can be accomplished by adding an extra tap to thefeedforward filter of the equalizer.

In another aspect of the invention, the simulated receiver uses a verylong equalizer that approximates the performance of an infinite lengthequalizer. This is useful to establish a performance bound on practicallength equalizers. In yet another aspect of the invention, a muchshorter equalizer is used in the simulated receiver. The shorterequalizer more closely approximates the performance of a practicallyrealizable receiver, but it requires additional optimization ofequalizer delay and sampling phase, examples of which are also describedherein.

Other aspects of the invention include methods and devices correspondingto the techniques described above.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention has other advantages and features which will be morereadily apparent from the following detailed description of theinvention and the appended claims, when taken in conjunction with theaccompanying drawings, in which:

FIG. 1 (prior art) is a block diagram of a 10 G optical fibercommunications link.

FIG. 2 is a block diagram of a communications link suitable for use withthe invention.

FIGS. 3 and 4 are block diagrams of a test system according to theinvention.

FIG. 5 is a block diagram of a model of a test system according to theinvention.

FIG. 6 is block diagram of a feedforward filter with a tap to adjust aconstant offset.

FIGS. 7A and 7B are block diagrams of mathematically equivalent modelsfor the sampler output.

FIG. 8 is a time plot illustrating calculation of optical modulationamplitude from a square wave that is time-aligned with an expectedtransition time.

FIGS. 9A and 9B are time plots illustrating miscalculation of opticalmodulation amplitude from square waves that are temporally misalignedwith an expected transition time.

FIG. 9C is a time plot illustrating calculation of optical modulationamplitude from a square wave that is temporally aligned with an expectedtransition time.

FIG. 10 is a flow chart of a method according to the invention fortemporally aligning a square wave synthesized from a captured waveform.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

FIG. 2 is a block diagram of a communications link 300 suitable for usewith the current invention. The link 300 includes a transmitter 305coupled through a channel 310 to a receiver 320. The overall system 300may suffer from various impairments in transmitting data from thetransmitter 305 to the receiver 320. The transmitters 305, channels 310and receivers 320 may be various types, depending on the endapplication. Microwave, RF, cable, and optical fiber are some examplesof different media for the channel 310. Different types of modulationmay be used by the transmitter 305. Some examples are on-off keying,QAM, PSK, and OFDM. Similarly, the receiver 315 can also take thevarious forms corresponding to these different modulation formats.

FIGS. 3 and 4 are block diagrams of a test system according to thecurrent invention. The purpose of the test system is to test thetransmitter, for example for compliance with a specific standard. InFIG. 3, a data test pattern 450 is applied to the transmitter under test405. The transmitter 405 produces an output 460, which is captured by awaveform capture device 465, for example a sampling oscilloscope. Thecaptured waveform may be stored in some medium 467, for subsequentretrieval.

In FIG. 4, a software simulator 470 accesses the stored waveform 462.The simulator 470 is designed to simulate the effects of a referencechannel 472 and/or reference receiver 474. The software 470 processesthe waveform 462 to simulate propagation of the waveform through thereference channel 472 and/or the reference receiver 474. A performancemetric 476 is calculated based on the simulated propagation. The datatest pattern 450 may be used as an input to the software simulator 470,for example in modeling the reference receiver/channel or in calculatingthe performance metric. By using this software simulator-based approach,the performance of the transmitter under test 405 can be predictedwithout requiring the use of hardware reference components.

Described below is a specific embodiment using the 10 GBASE-LRM systemas a specific example. The following example describes a testmethodology for measuring a performance metric known as the TransmitterWaveform and Dispersion Penalty (TWDP). TWDP is defined as

TWDP=SNR_(REF)−SNR_(TX)(dBo)   (1)

where SNR_(REF) is a reference signal to noise ratio (SNR), expressed inoptical dB (dBo), and SNR_(TX) is the equivalent SNR in dBo for thetransmitter under test at the slicer input of a reference decisionfeedback equalizer (DFE) receiver for the measured waveform afterpropagation through a reference optical fiber channel. The transmitterunder test is compliant with the standard if its computed TWDP is lessthan the upper limit specified by the IEEE 802.3aq standard.

In this test, the reference SNR, SNR_(REF), has a specific numericvalue, and it is related to OMA_(RCV), T, and N₀ through

SNR_(REF)=OMA_(RCV)√{square root over (T/(2N ₀))}

SNR_(REF)(dBo)=10 log₁₀(OMA_(RCV)√{square root over (T/(2N ₀)))}  (2)

where OMA_(RCV) is the OMA of the signal at the input to the referencechannel, N₀ is the one-sided power spectral density of the additivewhite Gaussian noise (AWGN) assumed for the reference DFE receiver, andT is the bit duration, also called the unit interval.

Eqn. (2) gives the SNR that would be realized by an ideal matched filterreceiver if the received waveform were an ideal rectangular NRZ waveformwith amplitude OMA_(RCV) and the receiver had AWGN of spectral densityN₀. This ideal received waveform would result, for example, if an idealNRZ transmitter transmitted an optical signal of amplitude OMA_(RCV)over a distortion-free, lossless optical fiber.

The specific value of SNR_(REF) is set equal to some margin in opticaldB above SNR_(RQD), where SNR_(RQD) is the SNR required for the idealmatched filter receiver to achieve a desired bit error rate (BER) for anideal NRZ received waveform. For the 10 GBASE-LRM example, the desiredBER is 10⁻¹², and SNR_(REF) is set 6.5 dBo above SNR_(RQD) to allow for6.5 dB of optical power penalties related to dispersion. The BER for asuch a matched filter receiver with an ideal NRZ input waveform is givenby

BER=Q(SNR_(RQD))   (3)

where Q( ) is the Gaussian error probability function

$\begin{matrix}{{Q(y)} = {\int_{y}^{\infty}{\frac{1}{\sqrt{2\pi}}^{- \frac{x^{2}}{2}}{x}}}} & (4)\end{matrix}$

For the 10 GBASE-LRM example, Eqns. (3) and (4) result in 10log₁₀(SNR_(RQD)) equal to 8.47 dBo to achieve a BER of 10⁻¹². Allowingfor 6.5 dB of optical power penalties related to dispersion results inan SNR_(REF) given by

SNR_(REF)=8.47+6.5=14.97 dBo   (5)

For a given OMA_(RCV), SNR_(REF) determines a particular N₀ through Eqn.(2). Without loss of generality, OMA_(RCV) is normalized to 1, and N₀ isset accordingly.

Turning now to the transmitter under test, in a preferred embodiment,the SNR_(TX) (and hence, also the TWDP) is calculated by MATLAB code asillustrated in FIG. 5. The data test pattern {x(n)} 550 driving thetransmitter under test 505 is a periodic PRBS9 or similar data pattern.The transmitter 505 outputs an optical waveform 560, which is convertedto an electrical signal by an optoelectronic converter 561, filtered bya hardware reference filter 562, and then digitized and captured 563. Inone embodiment, the optoelectronic conversion, hardware filtering, anddigitization and capture are all performed by a sampling oscilloscope565 with an appropriate optical head, though the functions mayalternatively be performed by separate devices. The oscilloscope is setto capture at least one complete cycle of the periodic signal 560 withat least seven samples per unit interval. The number of samples per unitinterval is not critical as long as the sampling rate is high enough tocapture the high frequency content of the signal out of the hardwarefilter 562 without aliasing. In a preferred embodiment, the hardwarefilter 562 has a fourth-order Bessel-Thomson response with 3 dBelectrical bandwidth of 7.5 GHz to filter the incoming waveform beforedigitization. The oscilloscope is set to average out noise in thewaveform, thus producing the digitized waveform 566 which approximates anoiseless filtered version of the transmitter output 560. In analternate embodiment, no hardware filter 562 is used. However, thisalternate embodiment may lead to variation in measurement from one testsetup to another depending of the inherent bandwidth of theoptoelectronic conversion and subsequent digitization.

The inputs to the software simulator 570 include the following:

-   -   The captured waveform 568, which, in this specific example, is        one complete cycle of the digitized waveform 566 with an integer        number of samples per bit period or unit interval—e.g., one        complete cycle of the waveform resulting from a periodic PRBS9        data test pattern. In some cases, the digitized waveform 566 is        pre-processed (e.g., re-sampled, truncated and/or aligned with        the data test pattern, as necessary) to produce the captured        waveform 568. In other cases, pre-processing may not be        necessary and the digitized waveform 566 can be used directly as        the captured waveform 568. In a preferred embodiment, the        captured waveform 568 has 16 samples per unit interval. The        number of samples per unit interval is not critical as long as        the sampling rate is high enough to represent the high frequency        content of the signal without aliasing.    -   One complete cycle of the data test pattern 550 used to generate        the captured waveform 568. The data test pattern 550 and the        captured waveform 568 are aligned (i.e., a rectangular pulse        train based on the data test pattern 550 is aligned with the        captured waveform 568 within one unit interval). The data test        pattern 550 has period N (e.g., 511 for PRBS9) and one cycle is        denoted {x(n)}, where 0≦n≦N-1.        The preprocessing steps of re-sampling, truncation to one        complete cycle, and alignment with the data test pattern are        shown performed in block 567 outside the simulator 570 and        outside the sampling oscilloscope 565. Alternatively, some or        all of these steps may be performed as part of the simulator        570, or inside the sampling oscilloscope 565, or as part of the        digitization and capture 563.

In a preferred embodiment, the captured waveform 568 is processed by thesoftware simulator 570, as follows:

1. The captured waveform 568 is normalized 571. The OMA and baseline(steady state zero level) of the waveform 568 are either measuredseparately and entered as inputs to the algorithm, or alternatively theyare estimated from the captured waveform, as described more fully below.The zero-level (baseline) is subtracted from the waveform 568 and thewaveform is scaled such that the resulting OMA is 1. Thus, thenormalized waveform has a baseline of 0 and an OMA of 1. N₀ is set suchthat SNR_(REF) is 14.97 dBo, as described above.

2. Propagation of the normalized waveform through an optical fiberreference channel 572 is simulated. In a preferred embodiment, thereference channel is represented by an impulse response. An alternativeembodiment can represent the reference channel by a frequency response,or by a lookup table, or by other mappings of input waveform to outputwaveform. In a preferred embodiment, propagation is accomplished in thefrequency domain. This is accomplished by performing Fourier transformson the normalized waveform and fiber impulse response, multiplying thetransforms together, and then performing an inverse Fourier transformback to the time domain. An alternative embodiment accomplishespropagation in the time domain by convolving the normalized waveformwith the fiber impulse response.

Note that the OMA can alternatively be normalized after propagationthrough the fiber, depending on the performance metric desired. As longas the fiber model is normalized to preserve OMA, normalizing before orafter propagation will give equivalent results.

3. The output from the reference channel 572 is passed through anantialiasing filter 573. A fourth-order Butterworth filter of bandwidth7.5 GHz is used in this example. An alternative embodiment would omitthe antialiasing filter 573 but would instead filter the white Gaussiannoise 577 before adding it to the reference fiber output 576.

4. The output signal from the antialiasing filter 573 is sampled 574 atrate 2/T with sampling phase φ. An alternate embodiment using a matchedfilter and a rate 1/T sampler is described below. The sampler output 575is denoted y_(φ)(nT/2), which has a deterministic component and a randomcomponent given by

y _(φ)(nT/2)=r(n)+η(n).

The sequence {r(n)} is the deterministic sampled version of the filteredoutput of the reference fiber; it is periodic with period 2N. Thesequence {η(n)} is a discrete-time noise sequence which is obtained bypassing the AWGN through the anti-aliasing filter and sampling.

5. The sampled signal is processed by a fractionally-spaced MMSE-DFEreceiver 584 with N_(f) feedforward taps (at T/2 spacing) 585 and N_(b)feedback taps 586. For example, N_(f) can be set to 14 and N_(b) can beset to 5. The use of 14 feedforward taps and 5 feedback taps are examplesettings; the number of taps can be different depending on the desiredcapabilities of the reference DFE receiver. The feedforward filter isoptionally augmented with an extra (N_(f)+1^(th)) tap that adjusts tooptimize a constant offset, as shown in FIG. 6 and described more fullybelow. The feedforward and feedback tap coefficients are calculatedusing a least-squares approach that minimizes the mean-squared error atthe slicer 578 input for the given captured waveform, assuming the noisepower spectral density N₀ set according to Eqn. (2) as described above.

FIG. 5 shows the channel and equalizer model used for the least-squarescalculation. The reference DFE consists of a feedforward filter {W(0), .. . , W(N_(f)−1)}, an optional constant offset coefficient W(N_(f)) (asshown in FIG. 6), and a feedback filter {B(1), . . . , B(N_(b))}. Aconventional DFE would feedback decisions {{circumflex over (x)}(n)}. Inthis case, the decisions are assumed to be correct and the decided bitsare therefore replaced with the transmitted bits {x(n)}. While a finitelength equalizer is shown in this example, an alternate embodiment wouldinclude an infinite length equalizer, performance of which could becomputed analytically.

The feedback filter 586 is symbol spaced and strictly causal, feedingback the N_(b) bits prior to the current bit x(n). The input sequence onwhich the slicer makes decisions is denoted {z (n)}, where

$\begin{matrix}\begin{matrix}{{z(n)} = {{\sum\limits_{k = 0}^{N_{f} - 1}{{W(k)}{y_{\varphi}\left( {{n\; T} + {DT} - {{kT}/2}} \right)}}} + {W\left( N_{f} \right)} -}} \\{{\sum\limits_{k = 1}^{N_{b}}{{B(k)}{x\left( {n - k} \right)}}}} \\{= {{\sum\limits_{k = 0}^{N_{f} - 1}{{W(k)}\left( {{r\left( {{2n} + {2D} - k} \right)} + {\eta \left( {{2n} + {2D} - k} \right)}} \right)}} + {W\left( N_{f} \right)} -}} \\{{\sum\limits_{k = 1}^{N_{b}}{{B(k)}{x\left( {n - k} \right)}}}}\end{matrix} & (6)\end{matrix}$

Eqn. (6) includes the term W(N_(f)), which is the optional constantoffset coefficient for automated offset compensation; without offsetcompensation, W(N_(f)) would be absent (or set to 0). In Eqn. (6), D isan integer such that the number of anticausal T/2-spaced taps in thefeedforward filter is 2D. While the feedforward filter is modeled ashaving anticausal taps, a practical equalizer would delay decisions byDT so that the overall system is causal. DT is referred to as theequalizer delay. Sometimes the equalizer delay is expressed as thenumber of feedforward anticausal taps, 2D.

The least-squares solution for the feedforward and feedback filtersminimizes the quantity

${M\; S\; E} = {E\left\lbrack {\sum\limits_{n = 0}^{N - 1}\left( {{z(n)} - {x(n)}} \right)^{2}} \right\rbrack}$

where the expectation operator, E, refers to the random noise at theslicer input, which results from the additive white Gaussian noise 577,filtered through the antialiasing filter 573, sampled 574, and filteredby the feedforward filter 585. (Here MSE is actually N times themean-squared error when averaged over the input sequence {x(n)}.) Thealgorithm finds the optimal W vector, B vector, and equalizer delay DTsuch that the mean-squared error is minimized in an efficient manner,for example as described in P. A. Voois, I. Lee, and J. M. Cioffi, “Theeffect of decision delay in finite-length decision feedbackequalization,” IEEE Transactions on Information Theory, vol. 53, pp.618-621, March 1996, which is incorporated herein by reference. Theoptimal sampling phase φ for the T/2 sampler is found by a brute-forcesearch over the 16 sampling phases in the unit interval. Note that for asufficiently long equalizer, that is, for N_(f) and N_(b) sufficientlylarge, the MSE is insensitive to sampling phase and equalizer delay. Inthat case, the number of anticausal taps 2D can be set to N_(f)/2, forexample, and optimization over D and φ can be omitted. Further detailsof the mean squared error (MSE) minimization are given below.

6. Once the MSE is minimized by finding the optimal sampling phase,equalizer delay, and tap coefficients, the bit-error rate is calculatedby the semi-analytic method. In an actual receiver, the slicer wouldmake binary decisions on its input based on whether the input is greaterthan or less than the slicer threshold. The probability of making anerror depends on the amplitude of the input signal, the statistics ofthe input noise, and the threshold of the slicer, which for this exampleis set to ½. The overall bit-error rate is calculated as follows:

-   -   a. The Gaussian noise variance at the input to the slicer 578 is        calculated.    -   b. For each bit in the data test pattern, the equalized input        signal to the slicer 578 is calculated and the probability of        error calculated based on the amplitude of the input signal, the        variance of the input noise, and the threshold value.    -   c. The probabilities of error are averaged over all bits in the        data test pattern to compute a total probability of error        BER_(TX).

7. The equivalent SNR in optical dB is derived from the BER_(TX)according to

SNR_(TX)=10 log₁₀(Q⁻¹(BER_(TX))) (dBo)   (7)

8. The TWDP for the simulated fiber reference channel 572 is equal tothe difference (in optical dB) between SNR_(REF) and the equivalentSNR_(TX), as shown in Eqn. (1) above.

There are additional alternate embodiments that are related to theexample algorithm described above. For example, the algorithm describedabove minimizes MSE at the slicer input, since that is how mostpractical equalizers are implemented. Alternatively, the referencereceiver could minimize bit error rate (BER) or some other performancemetric to compute the transmitter penalty. Minimization of BER would notnecessarily result in the same penalty as that computed by minimizingmean squared error.

As another example, the TWDP algorithm described above simulatespropagation through a single optical fiber reference channel.Alternatively, a number of optical fiber reference channels 572 may besimulated, each corresponding to a defined stressor that has beendefined to represent a certain level of dispersion and a certaincategory of fiber response. The normalized waveform is passed thougheach of these simulated reference channels 572 and the referenceequalizing receiver to compute a “channel-specific” TWDP for eachchannel. Transmitter compliance can be defined in terms of thecollection of channel-specific TWDP values. For example, compliance canrequire that a certain percentage of the channel-specific TWDP values bebelow some threshold. In particular, one could require that all of thechannel-specific TWDP values be below some limit. This is equivalent torequiring that the maximum of the channel-specific TWDP values be belowsome limit. Or, a different limit could be defined for differentreference channels or classes of reference channels. It will be readilyapparent that many different compliance criteria can be defined, all ofwhich should be considered within the scope of this invention.

Appendices A and B give examples of MATLAB code that illustrate variousaspects of the algorithm described above. The following derivationdetails relate the above description to the code presented in theappendices.

As described above, the TWDP algorithm requires estimates of the opticalmodulation amplitude (OMA) and baseline (steady state zero level) of thewaveform to normalize 571 the waveform. These values can be measuredexternally and input into the algorithm, which saves computation time,or the algorithm can estimate these quantities based on the inputwaveform. However, for distorted waveforms, such as those that mightresult from the use of lower cost optical transmitters, these values canbe difficult to measure accurately using the traditional method ofdisplaying a waveform on an oscilloscope. It is therefore advantageousfor the algorithm to automatically compute estimates of the OMA andbaseline from the input waveform 568 rather than relying on externallyprovided measurements.

One method of computing OMA and baseline, referred to as Method 1, is asfollows. In an optical communications system, the optical power outputby a laser is commonly modulated in a binary fashion to send data overan optical fiber. Nominally, the optical power is high for the durationof a bit period to send a logical “1”, and low to send a logical “0”.The difference in optical power between the nominal high level and thenominal low level is the OMA, and we refer to the nominal low level asthe baseline (or sometimes, the steady-state zero level, or simply thezero level). This type of modulation is commonly referred to ason-off-keying, where “on” means high laser power and “off” means lowlaser power. This is also referred to as nonreturn-to-zero (NRZ)signaling when the output power stays at nominally the same level for anentire bit period. In actuality, the level is not constant for theentire bit period due to finite rise and fall times of the laser output.Also, the laser output experiences ringing when it transitions from highpower to low power and from low power to high power. This ringing isreferred to as relaxation oscillation. Nonetheless, a common feature ofNRZ modulation is that a long string of zeros or a long string of oneswill each result in a modulated output signal that tends to a constantsteady state value.

A simple linear (strictly, affine) model of the laser output power y(t)is given by

$\begin{matrix}{{y_{lin}(t)} = {b + {O\; M\; A{\sum\limits_{n = {- \infty}}^{\infty}{{x(n)}{p\left( {t - {nT}} \right)}}}}}} & (8)\end{matrix}$

where b is a constant baseline value, x_(n)ε{0,1}, p(t) is the pulseshape output by the laser, and T is the bit period. OMA can be definedas the difference between the steady state “on” or “high” power of thelaser and the steady state “off” or “low” power of the laser. For thisdefinition of OMA to be consistent with Eqn. (8), a train of pulses p(t)must sum to 1, or

$\begin{matrix}{{\sum\limits_{n = {- \infty}}^{\infty}{p\left( {t - {nT}} \right)}} = 1} & (9)\end{matrix}$

If p(t) is the rectangular function Π_(T)(t), defined by

$\begin{matrix}{{\Pi_{T}(t)} = \left\{ \begin{matrix}1 & {{t} < {T/2}} \\{1/2} & {{t} = {T/2}} \\0 & {{t} > {T/2}}\end{matrix} \right.} & (10)\end{matrix}$

then y_(lin)(t) is ideal NRZ modulation.

A generalized NRZ modulation can be defined by allowing p(t) to be afiltered version of the rectangular pulse shape. Let h(t) be the impulseresponse of a filter normalized such that

$\begin{matrix}{{\int_{- \infty}^{\infty}{{h(t)}\ {t}}} = 1} & (11)\end{matrix}$

Let p(t)=Π_(T)(t)*h(t), where * denotes convolution. Then Eqn. (9) stillholds, as can be seen by

$\begin{matrix}\begin{matrix}{{\sum\limits_{n}{p\left( {t - {nT}} \right)}} = {{h(t)}*{\sum\limits_{n}{\Pi_{T}\left( {t - {nT}} \right)}}}} \\{= {{h(t)}*1}} \\{= {\int_{- \infty}^{\infty}{{h(t)}\ {t}}}} \\{= 1}\end{matrix} & (12)\end{matrix}$

Also,

∫_(−∞)^(∞)p(t) t = ∫_(−∞)^(∞)h(t) t∫_(−∞)^(∞)Π_(T)(t)t = T.

This is easily shown in the frequency domain.

Let q(t)=OMA p(t). Then Eqn. (8) is rewritten as

$\begin{matrix}{{y_{lin}(t)} = {b + {\sum\limits_{n = {- \infty}}^{\infty}{{x(n)}{q\left( {t - {nT}} \right)}}}}} & (13)\end{matrix}$

Now

$\begin{matrix}{{\int_{- \infty}^{\infty}{{q(t)}\ {t}}} = {{OMA} \cdot T}} & (14)\end{matrix}$

Hence, Eqn. (14) can be used to estimate OMA by finding q(t),integrating, and dividing by T. The pulse q(t) is referred to as theestimated pulse response of the transmitter or channel.

The problem of estimating OMA can now be cast as finding the b and q(t)in Eqn. (13) such that y_(lin)(t) approximates y(t). This problem issolved by finding the least-squares estimate that minimizes the meansquared error between y(t) and y_(lin)(t). This solution can beapproximated by transitioning to the discrete time domain and usinglinear algebra techniques.

Let t=mT+kΔ, where 0≦k≦K−1 and Δ=T/K Eqn. (13) can be generalizedslightly by allowing b to be a periodic function of t with period T.From Eqn. (13),

$\begin{matrix}\begin{matrix}{{y_{lin}\left( {{m\; T} + {k\; \Delta}} \right)} = {{b\left( {k\; \Delta} \right)} + {\sum\limits_{n = {- \infty}}^{\infty}{{x(n)}{q\left( {{\left( {m - n} \right)T} + {k\; \Delta}} \right)}}}}} \\{= {{b\left( {k\; \Delta} \right)} + {\sum\limits_{n = {- \infty}}^{\infty}{{x\left( {m - n} \right)}{q\left( {{nT} + {k\; \Delta}} \right)}}}}}\end{matrix} & (15)\end{matrix}$

Assume that the pulse response is of finite duration such that q(t)=0for t outside the interval [−AT,(M+1)T). Here A is the anticipation andM is the memory of the channel in bit periods.

Let y′_(m)≡[y(mT) y(mT+Δ) . . . y(mT+(K−1)Δ], with similar definitionsfor q_(m) and y_(lin,m). The notation v′ denotes the transpose of v. Letb′≡[b(0) . . . b((k−1)Δ)] and x′_(m)≡[x(m+A) . . . x(m) . . . x(m−M)].Eqn. (15) can be compactly expressed for all K sampling phases by

$\begin{matrix}\begin{matrix}{y_{{lin},m} = {b + {q_{- A}{x\left( {m + A} \right)}} + \ldots + {q_{M}{x\left( {m - M} \right)}}}} \\{= {b + {Qx}_{m}}}\end{matrix} & (16)\end{matrix}$

where Q is a K by A+M+1 matrix with columns q_(−A) through q_(M). Notethat each column of Q consists of K successive samples of thecontinuous-time pulse response q(t).

Now find the b and Q that minimize E[|ε|²], where

ε≡y−y_(lin,m)   (17)

|ε|² is the magnitude squared of the vector ε, and E [□] indicatesexpected value. Treat x and y as random vectors and drop the msubscript. The orthogonality principle for least squares estimationimplies that

E[(y−(b+Q*x))x′]=0   (18)

where Q* is the optimal Q for a given b. The least squares solutiongives

Q*=ρR _(x) ⁻¹ −bE[x′]R _(x) ⁻¹   (19)

where

R_(x)≡E[xx′]  (20)

and

ρ≡E[yx′]  (21)

Using Eqn. (19) in the expression for E[|ε²|] and minimizing withrespect to b gives the optimal b′as

$\begin{matrix}{b^{*} = \frac{\overset{\_}{y} - {\rho \; R_{x}^{- 1}\overset{\_}{x}}}{1 - {{\overset{\_}{x}}^{\prime}R_{x}^{- 1}\overset{\_}{x}}}} & (22)\end{matrix}$

where an overbar indicates expected value. Substituting b* for b in Eqn.(19) gives the optimal Q for a periodic baseline that is allowed to varywith sampling phase. To find the optimal Q for a baseline that isconstant, set b′=b[1 1 . . . 1] for a scalar b. Using this b and Q* fromEqn. (19) in the expression for E[|ε²|] and minimizing with respect to bleads to the optimal b*

$\begin{matrix}{b^{*} = {{1/K}{\sum\limits_{i = 1}^{K}b_{i}^{*}}}} & (23)\end{matrix}$

where b*_(i) is the i^(th) element of b*. Using b′=b*[1 1 . . . 1] inEqn. (19) gives the optimal Q for the constant baseline b*.

Using whichever optimal Q* is preferred (constant baseline or periodicbaseline), the samples of the continuous-time pulse responsecorresponding to this minimum mean squared error solution are given by

q((−A+m)T+kΔ)=Q* _(k+1,m+1) 0≦k≦K−1; 0≦m≦A+M   (24)

The integral in Eqn. (14) is approximated by

$\begin{matrix}{{\int_{- \infty}^{\infty}{{q(t)}\ {t}}} \approx {\Delta {\sum\limits_{k = 1}^{K}{\sum\limits_{m = 1}^{A + M + 1}Q_{k,m}^{*}}}}} & (25)\end{matrix}$

With T=KΔ, Eqns. (14) and (25) give

$\begin{matrix}{{OMA} \approx {\frac{1}{K}{\sum\limits_{k = 1}^{K}{\sum\limits_{m = 1}^{A + M + 1}Q_{k,m}^{*}}}}} & (26)\end{matrix}$

It turns out that the estimate of OMA is identical whether the Q* basedon constant baseline or periodic baseline is used in Eqn. (26). For theTWDP computation, the constant baseline is used as the estimatedbaseline.

The above results are now translated into MATLAB code to compute OMA forcaptured experimental data. Let y_(m) be the captured samples as definedabove, where each vector y_(m) consists of K samples of a single bitperiod of the sampled waveform for a repetitive data test pattern oflength N. Let x_(m) also be defined as above, where x_(m) is a vectorconsisting of A data bits preceding the “current” bit, the current bitx(m), and M data bits after the current bit. Form the matrices Y and Xwith columns y_(m) and x_(m) respectively, where m varies from 0 to N-1.With expectations taken over the N data vectors x_(m), each assumedequally likely, the expected values required for the least squares fitand OMA estimation above are given by

$\begin{matrix}{{R_{x} = {\frac{1}{N}X\; X^{\prime}}}{\rho = {\frac{1}{N}Y\; X^{\prime}}}{\overset{\_}{x} = {X\; {1/N}}}{\overset{\_}{y} = {Y\; {1/N}}}} & (27)\end{matrix}$

where 1 is an all-ones column vector of length N.

Eqn. (16) can be written compactly for all N bit periods as

$\begin{matrix}{Y_{lin} = {\left\lbrack {Q\mspace{14mu} b} \right\rbrack \begin{bmatrix}X \\1^{\prime}\end{bmatrix}}} & (28)\end{matrix}$

The optimal Q* and b* can be found using a pseudo-inverse and thecaptured waveform matrix Y:

$\begin{matrix}{\left\lbrack {Q^{*}\mspace{14mu} b^{*}} \right\rbrack = {{Y\left\lbrack {X^{\prime}\mspace{14mu} 1} \right\rbrack}\left\lbrack {\begin{bmatrix}X \\1^{\prime}\end{bmatrix}\left\lbrack {X^{\prime}\mspace{14mu} 1} \right\rbrack} \right\rbrack}^{- 1}} & (29)\end{matrix}$

With some manipulation and making use of the relations in Eqn. (27), itcan be shown that Eqn. (29) gives the same result as Eqn. (19) and Eqn.(22) for the periodic bias case. OMA and baseline can then be estimatedusing Eqns. (29), (23), and (26).

Eqns. (29), (23), and (26) are the basis of the following example MATLABcode, which can be used in the implementation of TWDP calculation toestimate OMA (MeasuredOMA in the code) and the baseline (SteadyZeroPowerin the code) from the input waveform 568. The variables ant and mem arethe anticipation and memory, respectively, in bit periods. Theseparameters can be varied to fit the memory span of the intersymbolinterference of the input waveform. The accuracy of the OMA estimate isdependent on the choice of ant and mem. If ant and mem are set too low,the entire pulse response will not be contained within the span. On theother hand, if they are chosen very large, the 2^(A+M+1) possible binaryvectors of length A+M+1 will be sparsely represented in the data testpattern if the sequence is of length significantly less than 2^(A+M+1).In that case, the expected values computed may be significantlydifferent than the expectations taken over truly random data.

  ant=4; mem=5; X=zeros(ant+mem+1,PRBSlength);Y=zeros(OverSampleRate,PRBSlength); X(1,:)=[XmitData(1+ant:end)′XmitData(1:ant)′]; for ind=2:ant+mem+1  X(ind,:)=[X(ind−1,end)X(ind−1,1:end−1)]; end X=[X;ones(1,PRBSlength)]; forind=1:OverSampleRate Y(ind,:)=yout0([0:PRBSlength−1]*OverSampleRate+ind)′; endQmat=Y*X′*(X*X′){circumflex over ( )}(−1);SteadyZeroPower=mean(Qmat(:,ant+mem+2)); Qmat=Qmat(:,1:ant+mem+1);MeasuredOMA=sum(Qmat(:))/OverSampleRate;

Method 1 for OMA and baseline estimation used a definition of OMA as thedifference between the steady state “on” or “high” power of the laserand the steady state “off” or “low” power of the laser. However, IEEE802.3ae Clause 52.9.5 specifies a method of measuring OMA based on atransmitted square wave pattern that can result in an OMA measurementdifferent from this definition. The reason is that the laser output maynot reach its steady state “on” value or steady state “off” value duringthe half period of the square wave that the laser is respectively “on”or “off”. The method described next, referred to as Method 2, results inan OMA value consistent with the measurement technique specified in IEEE802.3ae Clause 52.9.5.

Like the method presented previously, Method 2 performs a linear fitthat finds the pulse response of the modulation that gives the best fitto the captured waveform. Instead of integrating that pulse response toestimate the OMA, Method 2 uses that pulse response to synthesize whatthe optical waveform would have been if a square wave of the appropriatefrequency had been transmitted. The amplitude of the optical signalduring measurement windows is then measured and used to determine OMA.For example, the amplitude of the optical signal at the center of the“on” or “high” portion of the square wave is measured by averaging theamplitude over a measurement window positioned at the center 20% of the“on” portion, and the amplitude of the optical signal at the center ofthe “off” or “low” portion of the square wave is similarly measured overa measurement windows positioned at the center 20% of the “off” portion.The OMA estimate is then the difference between the “on” amplitude andthe “off” amplitude, and the baseline estimate is the measurement of the“off” amplitude.

IEEE 802.3ae Clause 52.9.1.2 specifies that the data test pattern forthe square wave may be “a pattern consisting of four to elevenconsecutive ones followed by an equal run of zeros”. A preferredembodiment uses a square wave comprising a repetitive pattern with eightconsecutive ones followed by eight consecutive zeros, though otherperiod lengths (for example, ten ones and ten zeros) could be usedinstead. An example implementation of Method 2 for OMA and baselineestimation is included in the MATLAB code in Appendix B.

FIG. 8 shows an example calculation of OMA using Method 2. In FIG. 8,the synthesized square waveform 810 includes eight consecutive zerosfollowed by eight consecutive ones. The OMA of the synthesized waveform810 is calculated as the difference between the average value of thesynthesized square waveform 810 captured during measurement window 820and the average value of the synthesized square waveform 810 capturedduring measurement window 830. Measurement windows 820, 830 are aligned,respectively, with the central regions of the consecutive zeros and theconsecutive ones. For example, measurement window 820 is aligned withthe center of the eight consecutive zeros and measurement window 830 isaligned with the center of the eight consecutive ones. In oneembodiment, the temporal position of the measurement windows 820,830 isspecified assuming an ideal square wave transitioning from zero to oneat a specified time. In the example shown in FIG. 8, the measurementwindows 820, 830 are temporally positioned based on an expectedtransition time from zero to one occurring at 8 unit intervals (UI).However, in other embodiments, the time when the ideal waveformtransitions between different values determines the temporal position ofthe measurement windows 820, 830. Because waveforms associated with thedirect output of the transmitter generally exhibit fast edge rates, thevalues of these waveforms typically are settled within the measurementwindows 820, 830, as shown in FIG. 8. While FIG. 8 shows ideal temporalalignment of the synthesized square waveform 810, some temporalmisalignment of the synthesized waveform 810 can be tolerated in certainimplementations. For example, due to the fast settling times, OMA may beaccurately calculated if the synthesized waveform 810 is aligned within1 unit interval (UI) of the expected transition time assumed in temporalpositioning of the measurement windows 820, 830.

Because the synthesized waveform 810 shown in FIG. 8 has stable values,it is relatively insensitive to misalignment between the synthesizedwaveform and the measurement windows. However, other synthesizedwaveforms may be very sensitive to timing misalignment. For example, OMAmay be calculated after propagation through a channel having dispersiveeffects. The dispersive effects of the channel can result in a widerpulse response. The corresponding synthesized waveform can have moregradually sloping edges, possibly consuming much or all of the run ofconsecutive zeros and/or the run of consecutive ones. In this case, atemporal misalignment between the synthesized waveform and measurementwindows reduces accuracy of OMA calculation and the resulting TWDPcalculation.

FIGS. 9A and 9B illustrate this effect. In these examples, thesynthesized waveforms 910A, 910B are unsettled during the measurementwindows 920, 930. For example, in FIG. 9A, the synthesized waveform 910Ahas not settled on a zero value during measurement window 920 and alsohas not settled on a one value during measurement window 930. Rather,the waveform is still sloping within the measurement windows.Furthermore, the synthesized square waveform 910A is temporallymisaligned with an expected transition time which temporally positionsthe measurement windows 920, 930. In the examples of FIGS. 8A and 9B,the expected transition time is 8 unit intervals (UI), so themeasurement windows 920, 930 in FIG. 9A are centered at 4 and 12 UI.FIG. 9C shows an example of a temporally aligned waveform 910C having azero crossing at 8 UI. However, waveform 910A has a zero crossing atapproximately 7, exhibiting an approximately 1 UI misalignment. This 1UI misalignment also causes synthesized waveform 910A to include anadditional zero crossing at approximately 15 UI. The temporalmisalignment and unsettled waveform mean that the values within themeasurement windows 920, 930 are shifted up or down, resulting in errorsin the values measured within the measurement windows 920, 930. Theseerrors in turn can result in an error in the calculated OMA.

For example, the temporally misaligned waveform 910A in FIG. 9A producesa calculated OMA of 1.8801 rather than the true OMA of 1.7373 producedby the temporally aligned waveform 910C shown in FIG. 9C. Similarly, inFIG. 9B, a temporal misalignment of the synthesized waveform 910B duringthe measurement windows 920, 930 produces an OMA of 1.545 rather thanthe true OMA of 1.7373 calculated from the temporally aligned waveform910C.

Because measurement windows are temporally positioned based on expectedtiming of a synthesized waveform, such as an expected transition timewhen the synthesized waveform transitions between logic levels, temporalmisalignment between an actual synthesized waveform and the expectedtiming of the synthesized waveform alters the values obtained within themeasurement windows which decreases OMA calculation accuracy. Forexample, if the actual synthesized waveform transitions between logiclevels at a time different than the expected transition time, the actualsynthesized waveform is temporally misaligned with the measurementwindows. To improve OMA calculation accuracy and improve TWDPcalculation, temporal misalignment between the synthesized waveform andits expected timing, such as the expected transition time, is correctedby time-shifting the synthesized received waveform relative to themeasurement windows. FIG. 10 is a flow diagram illustrating a method fortemporally aligning the synthesized received waveform and themeasurement windows. The temporal alignment method described below maybe used to calculate the OMA of a waveform captured before transmissionthrough a physical channel or of a waveform captured at the output of aphysical channel having dispersive effects.

A waveform is initially captured 1005 and a pulse response is derived1020 from the captured waveform. For example, a transmitter generates awaveform corresponding to an applied data test pattern, such as aperiodic PRBS 9 signal or similar data pattern and the generatedwaveform is captured 1005 either directly from a transmitter or afterpropagation through a channel. The pulse response is then derived 1020from the captured waveform. In other implementations, the capturedwaveform is propagated through a simulated channel and the pulseresponse is derived 1020 from a waveform at the output of the simulatedchannel.

The derived pulse response is then applied to an ideal square wave tosynthesize 1030 a waveform for OMA calculation. In one example, thesquare wave is a sequence of eight consecutive zeros followed by eightconsecutive ones. The synthesized waveform may be temporally misalignedwith the measurement windows, which are temporally positioned based onan expected transition time between logic levels. For example, themeasurement windows shown in FIGS. 8, 9A, 9B and 9C are centered at 4unit intervals (UI) and 12 UI based on an expected transition time of 8UI. If the synthesized waveform has a different transition time, such as9 UI, the synthesized waveform and measurement windows are misaligned.

To compensate for temporal misalignment, a specific transition between afirst logic level (such as a “zero”) and a second logic level (such as a“one”) is located 1040 in the synthesized waveform. In oneimplementation, the synthesized waveform is AC-coupled, so thetransition is initially approximated as the time when the waveform has avalue of zero, simplifying transition location 1040. In otherimplementations, the synthesized waveform may be similarly modified sothat transitions are initially approximated as the time when thewaveform has a specific value. To locate 1040 a specific transition,certain transitions may be ignored. For example, in FIGS. 9A-9B, thedesired specific transition is the central transition located near 8 UI.In this example, transitions located within a specified interval, suchas 2 UI, from the beginning and end of a waveform analysis region may bediscarded so that the central transition is properly identified. Thetemporal misalignment between the synthesized waveform and themeasurement windows is determined by the misalignment between thewaveform used to derive the pulse response and the reference datapattern. In the present example, the waveform used to derive 1020 thepulse response is aligned with the reference data pattern within 1 UI.This allows the timing of the synthesized waveform to match that of anideally-timed waveform, such as the synthesized waveform 810 shown inFIG. 8, within 1 UI. With this constraint, the desired centraltransition occurs with 1 UI of its ideal location, allowing transitionsoccurring within 2 UI of the beginning and end of the waveform analysisregion to be ignored. If the waveform is sampled, then the two samplesthat straddle the time the waveform has the specific value can beinterpolated to more precisely locate 1040 the time when the waveformtransitions between the first logic level and the second logic level.

The time when the identified 1040 transition between the first logiclevel and the second logic level occurred is compared 1050 to theexpected transition time which was used to temporally position themeasurement windows and the temporal offset between the transition inthe synthesized waveform and the expected transition time is determined.The synthesized waveform is temporally shifted 1060 by the determinedtemporal offset so the transition time of the synthesized waveform isaligned with the expected transition time, aligning the synthesizedwaveform with the measurement windows. OMA can then be estimated basedon the values obtained during the measurement windows and used tocalculate TWDP.

The following example MATLAB code can be used to implement the methoddescribed above to take in the temporally misaligned synthesizedwaveform (Y0 in the code) to produce an aligned synthesized waveform(Yin the code). The example MATLAB code can be used in an implementationof TWDP calculation to improve estimation of OMA. Appendix C providesadditional example MATLAB code implementing a TWDP calculation includingthe temporal alignment described above.

  function Y = AlignY(Y0,SqWvPer,OverSampleRate) % Aligns the crossingsof the xMA square waveform Y % to the first sample in the vector, etc.The first % sample is equivalent to time = 0. Y = Y0−mean(Y0); %AC-couple so crossings are at 0. % Look only for the approximatecrossing in the % middle by ignoring any within ~2 UI from the %beginning. Due to possible misalignment of the % captured waveform, thisis the only crossing that % is certain. x =find(sign(Y(2*OverSampleRate:end−1)) ~= . . .sign(Y(2*OverSampleRate+1:end)),1)+2*OverSampleRate−1; % Use linearinterpolation to find a more exact % crossing point. xinterp =interp1([Y(x),Y(x+1)],[x,x+1],0); % Calculate the difference from theideal crossing % time, then shift and re-sample to create the % alignedsquare waveform. Y = circshift(Y0,SqWvPer/2*OverSampleRate+1−x); Y =interp1(Y, (1:length(Y))′+xinterp−x,′spline′);

Alternatively, the MATLAB code presented below describes an alternateimplementation of the time alignment method shown in FIG. 10 to take inthe temporally misaligned synthesized waveform (Y0 in the code) andproduce an aligned synthesized waveform (Y in the code). This exampleMATLAB provides an alternate technique for initially approximating thetransition between low logic level and high logic level in thesynthesized waveform.

  Y = Y0−mean(Y0); % AC-couple so crossings are at 0. tmp =Y(1:end−1).*Y(2:end); % Crossings will be % where this vector is <=0.Look only for the % approximate crossing in the middle by ignoring any %within ~2 UI from the beginning. Due to possible % misalignment of thecaptured PRBS waveform, this % is the only crossing that is certain. x =find(tmp(2*OverSampleRate:end) <=0,1) + 2*OverSampleRate−1; % Use linearinterpolation to find a more exact % crossing point. xinterp =x−Y(x)/(Y(x+1)−Y(x)); % Calculate the difference from the ideal crossing% time then shift and re-sample to create the % aligned square waveform.xadj = xinterp−(SqWvPer/2*OverSampleRate+1); Y =circshift(Y0,−floor(xadj)); Y = interp1(Y,(1:length(Y))″+xadj−floor(xadj),′spline′);

Accordingly, temporally aligning the synthesized waveform andmeasurement windows allows for more accurate calculation of OMA,particularly when the synthesized waveform is subject to distortion.Improving OMA calculation consequently enables more accurate evaluationof TWDP to better determine transmitter performance. Although the methoddescribed above addresses temporal alignment for OMA calculation, theabove method may also be used to improve calculation of modulationamplitude in a variety of units, such as Voltage Modulation Amplitude(VMA) calculation or Current Modulation Amplitude (CMA). Hence, allowingthe temporal alignment method also improve performance analysis, such aswaveform and dispersion penalty (WDP), depending on modulationamplitude, at any point in a system, including a receiver.

Automated offset optimization is now discussed in additional detail.While the offset optimization technique described here is presented inthe context of a simulated reference receiver for a transmitter test, itcan also be used in a real equalizing receiver to correct offset in thereceived waveform. In one embodiment, the TWDP algorithm corrects forinaccuracies in the estimated value of the baseline (or steady statezero level). This automated offset method makes the reported TWDP valueindependent of estimation error in the baseline. In fact, it obviatesthe need to measure a baseline at all provided the ratio of baselinevalue to OMA is not unreasonably large. After OMA normalization andsubtraction of the estimated baseline, the resulting normalized waveformhas target signal values of “0” and “1”. The slicer threshold for BERdetermination is set at ½. An error in baseline estimation such that theresulting waveform is improperly offset and not centered about ½ willresult in a higher BER than if the offset were optimized. This offsetoptimization can be alternatively viewed as a threshold optimization.

The automated offset optimization can be accomplished by extending thevector of feedforward equalizer tap coefficients by 1. The additionaltap is not fed by a delayed version of the received signal, but isinstead fed by the constant 1. FIG. 6 shows an example feedforwardfilter with N_(f) regular taps and an extra (N_(f)+1)^(th) tap foroffset optimization. The MMSE tap coefficients are then calculated inthe normal fashion, and the coefficient of the extra tap essentiallygives the offset that minimizes mean squared error with the targetsignal values. An alternate method extends the length of the feedbacktap coefficient vector by one, feeding that extra feedback tap by theconstant 1, and computing the offset as part of the feedback tapcoefficient vector.

A derivation of the MSE minimization and offset optimization is nowprovided that relates the mathematical description above to the sampleMATLAB code in Appendices A and B. FIGS. 7A and 7B show mathematicallyequivalent models for generating the sampler output sequencey_(φ)(nT/2). The anti-aliasing filters 573 are identical, as are the T/2samplers 574 (identical sampling rate and phase). The model in FIG. 7Bshows explicitly the random noise component η(n) and the deterministicperiodic component r(n) that comprise the sampler output y_(φ)(nT/2).The actual values of the sequence {η(n)} are not required for MSEminimization; only the second order statistics are required. Thissuggest an alternate embodiment where the noise sequence η(n) isdescribed only in terms of its second order statistics, and AWGN 577 isnot necessarily assumed in the model.

Assume that the AWGN 577 is mean zero, in which case {η(n)} is also meanzero. Denote the autocorrelation sequence for {η(n)} by {c(k)}.

c(k)≡E[η(n)η(n−k)]=c(−k)

The derivation for the case without offset optimization (W(N_(f))=0) isdescribed first. From Eqn. (6), the slicer input for deciding bit n is

$\begin{matrix}{{z(n)} = {{\sum\limits_{k = 0}^{N_{f} - 1}{{W(k)}\left( {{r\left( {{2n} + {2D} - k} \right)} + {\eta \left( {{2n} + {2D} - k} \right)}} \right)}} - {\sum\limits_{k = 1}^{N_{b}}{{B(k)}{x\left( {n - k} \right)}}}}} & (30)\end{matrix}$

This can be conveniently described in matrix notation for all values ofn, 0≦n≦N-1, as

$\begin{matrix}{\mspace{79mu} {{z = {{\left( {R + N} \right)w} - {X\begin{bmatrix}0 \\b\end{bmatrix}}}}\mspace{76mu} {where}}} & (31) \\\begin{matrix}{\mspace{76mu} {{z \equiv \begin{bmatrix}{z(0)} \\{z(1)} \\\vdots \\{z\left( {N - 1} \right)}\end{bmatrix}};}} & {{w \equiv \begin{bmatrix}{W(0)} \\{W(1)} \\\vdots \\{W\left( {N_{f} - 1} \right)}\end{bmatrix}};} & {b \equiv \begin{bmatrix}{B(1)} \\{B(2)} \\\vdots \\{B\left( N_{b} \right)}\end{bmatrix}}\end{matrix} & \left( {32a} \right) \\{R \equiv \begin{bmatrix}{r\left( {2D} \right)} & {r\left( {{2D} - 1} \right)} & \ldots & {r\left( {{2D} - N_{f} + 1} \right)} \\{r\left( {{2D} + 2} \right)} & {r\left( {{2D} + 1} \right)} & \ldots & {r\left( {{2D} - N_{f} + 3} \right)} \\\vdots & \vdots & \; & \vdots \\{r\left( {{2D} + {2\left( {N - 1} \right)}} \right)} & {r\left( {{2D} + {2N} - 3} \right)} & \ldots & {r\begin{pmatrix}{{2D} + {2\left( {N - 1} \right)} -} \\{N_{f} + 1}\end{pmatrix}}\end{bmatrix}} & \left( {32b} \right) \\{N \equiv {\quad\begin{bmatrix}{\eta \left( {2D} \right)} & {\eta \left( {{2D} - 1} \right)} & \ldots & {\eta \left( {{2D} - N_{f} + 1} \right)} \\{\eta \left( {{2D} + 2} \right)} & {\eta \left( {{2D} + 1} \right)} & \ldots & {\eta \left( {{2D} - N_{f} + 3} \right)} \\\vdots & \vdots & \; & \vdots \\{\eta \left( {{2D} + {2\left( {N - 1} \right)}} \right)} & {\eta \left( {{2D} + {2N} - 3} \right)} & \ldots & {\eta \begin{pmatrix}{{2D} + {2\left( {N - 1} \right)} -} \\{N_{f} + 1}\end{pmatrix}}\end{bmatrix}}} & \left( {32c} \right) \\{\mspace{76mu} {X \equiv \begin{bmatrix}{x(0)} & {x\left( {- 1} \right)} & \ldots & {x\left( {- N_{b}} \right)} \\{x(1)} & {x(0)} & \ldots & {x\left( {{- N_{b}} + 1} \right)} \\\vdots & \vdots & \; & \vdots \\{x\left( {N - 1} \right)} & {x\left( {N - 2} \right)} & \ldots & {x\left( {N - N_{b} - 1} \right)}\end{bmatrix}}} & \left( {32d} \right)\end{matrix}$

For a given equalizer delay DT and sampling phase φ, minimizing MSErequires finding the vectors w and b that minimize the quantity

MSE=E[|z−x| ²]  (33)

where |□|² indicates the squared magnitude of the vector argument. Herex is the vector of N bits in the periodic sequence {x(n)}

$\begin{matrix}{x \equiv \begin{bmatrix}{x(0)} \\{x(1)} \\\vdots \\{x\left( {N - 1} \right)}\end{bmatrix}} & (34)\end{matrix}$

Using the fact that E[N] is the all-zeroes matrix, it follows that

E[|z−x| ² ]=w′(R′R+C _(η))w−2w′R′X{tilde over (b)}+{tilde over(b)}′X′X{tilde over (b)}  (35)

where v′ indicates the transpose of the vector or matrix v, and

$\begin{matrix}{{C_{\eta} \equiv {E\left\lbrack {N^{\prime}N} \right\rbrack}} = {N\begin{bmatrix}{c(0)} & {c(1)} & \ldots & {c\left( {N_{f} - 1} \right)} \\{c(1)} & {c(0)} & \ldots & {c\left( {N_{f} - 2} \right)} \\\vdots & \vdots & \; & \vdots \\{c\left( {N_{f} - 1} \right)} & {c\left( {N_{f} - 2} \right)} & \ldots & {c(0)}\end{bmatrix}}} & (36)\end{matrix}$

Minimizing the MSE in Eqn. (35) with respect to w for a given {tildeover (b)} gives

w _(opt)({tilde over (b)})=(R′R+C _(η))⁻¹ R′X{tilde over (b)}  (37)

Note that this result gives the optimal weight vector for a feedforwardonly equalizer by setting N_(b)=0. In that case {tilde over (b)}=1 andX=x. This can be used in an alternate embodiment that uses a linear(feedforward only) equalizer instead of a DFE in the reference receiver.

Substituting the expression for w_(opt)({tilde over (b)}) for w in Eqn.(35) gives

$\begin{matrix}\begin{matrix}{{E\left\lbrack {{z - x}}^{2} \right\rbrack} = {{\overset{\sim}{b}}^{\prime}\overset{\sim}{P}\overset{\sim}{b}}} \\{= {{\left\lbrack {1\mspace{14mu} b^{\prime}} \right\rbrack \begin{bmatrix}P_{00} & p^{\prime} \\p & P\end{bmatrix}}\begin{bmatrix}1 \\b\end{bmatrix}}} \\{= {P_{00} + {2p^{\prime}b} + {b^{\prime}P\; b}}}\end{matrix} & (38)\end{matrix}$

where

{tilde over (P)}≡X′X−X′R(R′R+C_(η))⁻¹R′X   (39)

Minimizing the resulting expression for MSE with respect to b gives

b _(MMSE) =−P ⁻¹ p   (40)

Collecting the results for the final answer, the minimizing feedbackvector b is

b _(MMSE) =P ⁻¹ p   (41 a)

where P is the square matrix in the lower right corner of {tilde over(P)} that excludes the first column and first row of {tilde over (P)}, pis the first column of {tilde over (P)} excluding the top element P₀₀,and {tilde over (P)} is given by

{tilde over (P)}≡X′X−X′R(R′R+C_(η))⁻¹R′X   (41b)

The minimizing feedforward vector w is given by

$\begin{matrix}{w_{MMSE} = {\left( {{R^{\prime}R} + C_{\eta}} \right)^{- 1}R^{\prime}{X\begin{bmatrix}1 \\b_{MMSE}\end{bmatrix}}}} & \left( {41c} \right)\end{matrix}$

The minimum MSE is given by

$\begin{matrix}\begin{matrix}{{MMSE} = {P_{00} + {2p^{\prime}b_{MMSE}} + {b_{MMSE}^{\prime}P\; b_{MMSE}}}} \\{= {P_{00} - {p^{\prime}P^{- 1}p}}}\end{matrix} & \left( {41d} \right)\end{matrix}$

The variance of the Gaussian noise at the input to the slicer, used tocompute the BER semianalytically, is given by

$\begin{matrix}{\sigma^{2} = {\frac{1}{N}w^{\prime}C_{\eta}w}} & (42)\end{matrix}$

Including the extra feedforward tap for automatic offset optimizationrequires a modification to the derivation above. Augment the feedforwardvector w with an extra tap to give the new feedforward vector {tildeover (w)}

$\begin{matrix}{\overset{\sim}{w} = \begin{bmatrix}w \\{W\left( N_{f} \right)}\end{bmatrix}} & (43)\end{matrix}$

Augment the matrix of noiseless channel outputs R with an additionalcolumn of all ones, denoted as the column vector 1, to give the newmatrix

{tilde over (R)}=[R 1]  (44)

Augment the matrix of noise samples N with an additional column of allzeros, denoted as the column vector 0, to give the new matrix

Ñ=[N 0]  (45)

Then

$\begin{matrix}{{{\overset{\sim}{C}}_{\eta} \equiv {E\left\lbrack {{\overset{\sim}{N}}^{\prime}\overset{\sim}{N}} \right\rbrack}} = \begin{bmatrix}C_{\eta} & 0 \\0^{\prime} & 0\end{bmatrix}} & (46)\end{matrix}$

The solution for the case with automated offset optimization is obtainedby replacing w with {tilde over (w)}, R with {tilde over (R)}, and C_(η)with {tilde over (C)}_(η) in Eqns. (41) and 42). Note that

$\begin{matrix}{{\left( {{{\overset{\sim}{R}}^{\prime}\overset{\sim}{R}} + {\overset{\sim}{C}}_{\eta}} \right)^{- 1} = \begin{bmatrix}{{R^{\prime}R} + C_{\eta}} & {R^{\prime}1} \\{1^{\prime}R} & N\end{bmatrix}^{- 1}};{{{\overset{\sim}{w}}^{\prime}{\overset{\sim}{C}}_{\eta}\overset{\sim}{w}} = {w^{\prime}C_{\eta}w}}} & (47)\end{matrix}$

which are used in the example MATLAB code in the appendices.

The alternate embodiment mentioned above that extends the feedbackfilter instead of the feedforward filter can be implemented, forexample, by extending the vector b by one element and appending a columnof all ones to the matrix X in Eqns. (41a) through (41d).

Appendix A includes sample MATLAB code that computes TWDP for areference DFE receiver with 100 T/2-spaced feedforward taps, 50 feedbacktaps, and automated offset optimization. The OMA and baseline of thesignal are entered manually instead of estimated automatically. Theequalizer delay is 50 feedforward taps (25 T), and the equalizer is longenough that optimization over equalizer delay and sampling phase are notrequired.

In the embodiment illustrated in Appendix A, the reference receiver 584uses a very long equalizer, closely approximating performance of anideal infinite length equalizer. The use of 100 feedforward taps 585 (atT/2 spacing), 50 feedback taps 586 (at T spacing), and an equalizerdelay of 50 feedforward taps (25 T) are representative values thatapproximate an infinite length equalizer for the optical fiber channelsof interest. Because the longer T/2 spaced feedforward equalizer isinsensitive to sample phase and exact value of equalizer delay, thisembodiment has the advantage that it does not require optimization overthese parameters. However, such a long equalizer may give overlyoptimistic results for certain channels and thus may not accuratelyreflect the equalizability of the waveform by a practical receiver thatuses a conventional DFE architecture. For example, because of the use ofpseudo-random sequences, certain non-linear distortions will appear asdelayed copies of the sequence added to the signal. If the delayed copyfalls within the span of the very long equalizer, the reference DFE willappear to correct this distortion. With truly random data (instead of aPRBS), the nonlinear distortion would not be equalizable by a classicDFE. Since the delayed copy is more likely to fall outside the span of ashorter equalizer, use of the shorter equalizer in the referencereceiver 584 is more likely to give a realistic indication of theequalizability of the signal in real receivers. Another example wherethe very long equalizer might give overly optimistic performancepredictions is when distortions in the transmitted waveform includeelectrical reflections that are delayed so much that they fall outsidethe span of a practical length equalizer but are effectively equalizedby the extremely long equalizer.

For these reasons, it is sometimes desirable to implement a much shorterDFE in the reference receiver, for example, an equalizer with 14feedforward taps and 5 feedback taps. The embodiment with a shortequalizer requires optimization of sampling phase and equalizer delay.This could be accomplished by setting the delay and phase to givenvalues, minimizing MSE over w and b as described above, and thenrepeating the process in a brute force manner until the optimum settingsfor equalizer delay and sampling phase are found. However, it issometimes possible to save computing time by performing an efficientoptimization over equalizer delay, depending on the equalizer length andthe search range for equalizer delay. One such method for efficientoptimization of equalizer delay is described in P. A. Voois, I. Lee, andJ. M. Cioffi, “The effect of decision delay in finite-length decisionfeedback equalization,” IEEE Transactions on Information Theory, vol.53, pp. 618-621, March 1996. This method is illustrated in the sampleMATLAB code in Appendix B.

In another alternative embodiment of the TWDP calculation, the receiverhas a matched filter, matched to the shape of the received pulse, whichallows the use of a T-spaced feedforward filter instead of a T/2-spacedfilter. For this alternative, the simulation algorithm estimates achannel pulse response based on the reference fiber output 576. Thelinear fit described above as Method 1 for OMA and baseline estimationcan be used, for example, to estimate the pulse response. The front endof the reference receiver then uses a matched filter matched to thispulse response. This matched filter replaces the antialiasing filter573. The signal is passed through the matched filter and is sampled atinteger multiples of T instead of T/2. However, this embodiment has thedisadvantage that the estimated pulse response may inaccurately modelthe true pulse shape due to nonlinearities in the transmitter and/orchannel. This may result in suboptimal equalization and couldunderestimate the link performance with a real receiver based on afractionally spaced equalizer.

Although the detailed description contains many specifics, these shouldnot be construed as limiting the scope of the invention but merely asillustrating different examples and aspects of the invention. It shouldbe appreciated that the scope of the invention includes otherembodiments not discussed in detail above. For example, a T-spacedequalizer without a matched filter could be used instead of a T/2-spacedequalizer if it is known that the simulated waveform entering thereference receiver will always have bandwidth less than 1/(2T). Asfurther examples, the specific number of taps used for feedforward andfeedback, and the anticipation and memory used in the linear fits shouldbe considered exemplary and not in any way limiting the scope of theinvention. Various other modifications, changes and variations whichwill be apparent to those skilled in the art may be made in thearrangement, operation and details of the method and apparatus of thepresent invention disclosed herein without departing from the spirit andscope of the invention as defined in the appended claims. Therefore, thescope of the invention should be determined by the appended claims andtheir legal equivalents.

APPENDIX A SAMPLE MATLAB CODE, EXAMPLE 1 This MATLAB code illustratesvarious aspects of the invention. This version implements a very longDFE that does not require optimization of equalizer delay or samplingphase. It includes code for automatic offset optimization and requiresexternal input of the OMA and steady-state zero level. %%%%%%% MATLAB(R) script to compute TWDP %%%%%%%%%%%%%%%%%%%%%% % %*** Includesmodification to automatically compute optimal offset making %*** errorsin SteadyZeroPower irrelevant. %% TWDP test inputs %% Units for alloptical power values must match. The %% values given below forTxDataFile, MeasuredWaveformFile, MeasuredOMA, %% and SteadyZeroPowerare examples and should be replaced by actual %% path\filenames andvalues for each waveform tested. %% Transmit data file: %% The transmitdata sequence is one of the TWDP test patterns. %% The file format isASCII with a single column of %% chronological ones and zeros with noheaders or footers. TxDataFile = ′prbs9_950.txt′; %% Measured waveform.The waveform consists of exactly K samples per unit %% interval T, whereK is the oversampling rate. The waveform must be %% circularly shiftedto align with the data sequence. The file format for %% the measuredwaveform is ASCII with a single column of chronological %% numericalsamples, in optical power, with no headers or footers.MeasuredWaveformFile = ′preproc-1207-01.txt′; %% OMA and steady-stateZERO power must be entered. MeasuredOMA = 6.6e−4; % Measured OMA, inoptical power SteadyZeroPower = 2.76e−4; % Measured steady-state logicZERO OverSampleRate = 16; % Oversampling rate %% Simulated fiberresponse, modeled as a set of ideal delta functions with %% specifiedamplitudes in optical power and delays in nanoseconds. %% FiberRespcontains delays in row 1, corresponding power coefficients in %% theother three rows. Rows 2, 3, and 4 each correspond to a differenct %%fiber channel model. For this example, use the power coefficients in %%row 4. The vector ′Delays′ contains the delays, and the vector 'PCoefs'%% contains the amplitudes. FiberResp = [. . .  0.000000 0.0727270.145455 0.218182  0.158 0.176 0.499 0.167  0.000 0.513 0.000 0.487 0.254 0.453 0.155 0.138]; Delays = FiberResp(1,:)′; PCoefs =FiberResp(4,:)′; %% Program constants %% SymbolPeriod = 1/10.3125; %Symbol period (ns) EFilterBW = 7.5; % Antialiasing filter bandwidth(GHz) %% For this example, use a very long equalizer, so optimization of%% equalizer delay and sampling phase are not required. Set theequalizer %% delay to half the length of the feedforward filter EqNf =100; % Number of feedforward filter taps EqNb = 50; % Number of feedbackfilter taps EqDel = ceil(EqNf/2); % Equalizer delay PAlloc = 6.5; %Allocated dispersion penalty (dBo) Q0 = 7.03; % BER = 10{circumflex over( )}(−12) %% Load input waveforms XmitData = load(TxDataFile); yout0 =load(MeasuredWaveformFile); PtrnLength = length(XmitData); TotLen =PtrnLength*OverSampleRate; Fgrid = [−TotLen/2:TotLen/2-1].′/(PtrnLength*SymbolPeriod); %% STEP 1 - NormalizeOMA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% yout0 = (yout0 - SteadyZeroPower)/MeasuredOMA; % Alternatively, we could just use yout0 =yout0/MeasuredOMA, since the % offset optimization below would correctfor failing to subtract off the % steady state zero level as long asSteadyZeroPower is not unreasonably % large compared to MeasuredOMA %%STEP 2 - Process waveform through simulated fiber channel. Fiber %%frequency response is normalized to 1 at DC (sum(PCoefs) =1). ExpArg =−j*2*pi*Fgrid; Hsys = exp(ExpArg * Delays′) * PCoefs; Hx =fftshift(Hsys/sum(PCoefs)); yout = real(ifft(fft(yout0).*Hx)); %% STEP3 - Process signal through front-end antialiasing filter %% Computefrequency response of front-end Butterworth filter [b,a] = butter(4,2*pi*EFilterBW,′s′); H_r = freqs(b,a,2*pi*Fgrid); %% Process signalthrough front-end filter yout = real(ifft(fft(yout) .* fftshift(H_r)));%% STEP 4 - Sample at rate 2/T %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% yout =yout(1:OverSampleRate/2:end); %% STEP 5 - Compute MMSE-DFE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% The MMSE-DFE filter coefficientscomputed below minimize mean-squared %% error at the slicer input. Thederivation follows from the fact that the %% slicer input over theperiod of the data sequence can be expressed as %% Z = (R+N)*W − X*[0B]′, where R and N are Toeplitz matrices constructed %% from the signaland noise components, respectively, at the sampled %% output of theantialiasing filter, W is the feedforward filter, X is a %% Toeplitzmatrix constructed from the input data sequence, and B is the %%feedback filter. The computed W and B minimize the mean square error %%between the input to the slicer and the transmitted sequence due to %%residual ISI and Gaussian noise. %% Compute the noise autocorrelationsequence at the output of the %% antialiasing filter and rate-2/Tsampler. N0 = SymbolPeriod/ (2 * Q0{circumflex over ( )}2 *10{circumflex over ( )}(2*PAlloc/10)); Snn = N0/2 *fftshift(abs(H_r).{circumflex over ( )}2) * 1/SymbolPeriod *OverSampleRate; Rnn = real(ifft(Snn)); Corr =Rnn(1:OverSampleRate/2:end); %% Construct a Toeplitz autocorrelationmatrix. C = toeplitz(Corr(1:EqNf)); %% Construct Toeplitz matrix frominput data sequence X = toeplitz(XmitData, [XmitData(1);XmitData(end:−1:end-EqNb+1)]); %% Construct Toeplitz matrix from signalat output of 2/T sampler. This %% sequence gets wrapped by equalizerdelay R = toeplitz(yout, [yout(1); yout(end:−1:end-EqNf+2)]); R =[R(EqDel+1:end,:); R(1:EgDel, :)]; R = R(1:2:end, :); %% Start ofmodifications to optimize offset ONE=ones(PtrnLength,1); %% Computeleast-squares solution for filter coefficients RINV =inv([R′*R+PtrnLength*C R′*ONE;ONE′*R ONE′*ONE]); R= [R ONE]; % Addall-ones column to compute optimal offset P = X′*X − X′*R*RINV*R′*X; P01= P(1,2:end); P11 = P(2:end,2:end); B = −inv(P11)*P01′; % Feedbackfilter W = RINV*R′*X*[1;B]; % Feedforward filter Z = R*W − X*[0;B]; %Input to slicer %% STEP 6 - Compute BER using semi-analytic method%%%%%%%%%%%%%%%%%%%%%% MseGaussian = W(1:end−1)′*C*W(1:end−1); %% End ofmodifications to optimize offset Ber =sum(0.5*erfc((abs(Z-0.5)/sqrt(MseGaussian))/sqrt(2)))/length(Z); %% STEP7 - Compute equivalent SNR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This %%function computes the inverse of the Gaussian error probability %%function. The built-in function erfcinv( ) is not sensitive enough for%% the low probability of error case. Q = inf; if Ber>10{circumflex over( )}(−12) Q = sqrt(2)*erfinv(1-2*Ber); elseif Ber>10{circumflex over( )}(−300) Q = 2.1143*(−1.0658−log10(Ber)).{circumflex over ( )}0.5024;end %% STEP 8 - Compute penalty %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%RefSNR = 10 * log10(Q0) + PAlloc; TWDP = RefSNR−10*log10(Q) % End ofprogram

APPENDIX B SAMPLE MATLAB CODE, EXAMPLE 2 This MATLAB code illustratesvarious aspects of the invention. This version implements a DFE with 14feedforward taps and 5 feedback taps. This code illustrates optimizationof equalizer delay and sampling phase, where the optimization ofequalizer delay is performed in an efficient manner. This versionincludes code that implements Method 2 for automatically estimating OMAand baseline from the input waveform and it includes automatic offsetoptimization This version computes TWDP for three different stressorchannels, using the maximum of those as the reported TWDP value. %%%%%%%MATLAB (R) script to compute TWDP %%%%%%%%%%%%%%%%%%%%%% %% TWDP inputs%% The values given below for TxDataFile and MeasuredWaveformFile are %%examples and should be replaced by actual path\filenames for each %%waveform tested. %% Transmit data file: The transmit data sequence isone of the TWDP test %% patterns. The file format is ASCII with a singlecolumn of chronological %% ones and zeros with no headers or footers.TxDataFile = ′prbs9_950.txt′; %% Measured waveform: The waveformconsists of exactly K samples per unit %% interval T, where K is theoversampling rate. The waveform must be %% circularly shifted to alignwith the data sequence. The file format for %% the measured waveform isASCII with a single column of chronological %% numerical samples, inoptical power, with no headers or footers. MeasuredWaveformFile =′preproc-1207-01.txt′; OverSampleRate = 16; % Oversampling rate, must beeven %% Simulated fiber responses, modeled as a set of ideal deltafunctions %% with specified amplitudes in optical power and delays innanoseconds. %% FiberResp contains delays in row 1, corresponding powercoefficients in %% the other three rows. Rows 2, 3, and 4 eachcorrespond to a differenct %% fiber channel model. For this example, aTWDP value is computed for %% each channel, and the TWDP reported is themaximum of the three TWDP %% values thus computed. The vector ‘Delays’contains the delays. FiberResp = [. . .  0.000000 0.072727 0.1454550.218182  0.158 0.176 0.499 0.167  0.000 0.513 0.000 0.487  0.254 0.4530.155 0.138+; Delays = FiberResp(1,:)′; %% Program constants %%SymbolPeriod = 1/10.3125; % Symbol period (ns) EqNf = 14; EqNb = 5; %Number of feedforward and feedback equalizer taps % Set search range forequalizer delay, specified in symbol periods. Lower % end of range isminimum channel delay. Upper end of range is the sum of % the lengths ofthe FFE and channel. Round up and add 5 to account for the %antialiasing filter. EqDelMin = floor(min(Delays)/SymbolPeriod);EqDelMax = ceil(EqNf/2 + max(Delays)/SymbolPeriod)+5; EqDelVec =[EqDelMin:EqDelMax]; PAlloc = 6.5; % Total allocated dispersion penalty(dBo) Q0 = 7.03; % BER = 10{circumflex over ( )}(−12) N0 =SymbolPeriod/2 / (Q0 * 10{circumflex over ( )}(PAlloc/10)){circumflexover ( )}2; EFilterBW = 7.5; % Antialiasing filter bandwidth (GHz) %%Load input waveform and data sequence, generate filter and other %%matrices yout0 = load(MeasuredWaveformFile); XmitData =load(TxDataFile); PtrnLength = length(XmitData); TotLen =PtrnLength*OverSampleRate; Fgrid =[−TotLen/2:TotLen/2-1].′/(PtrnLength*SymbolPeriod); % antialiasingfilter [b,a] = butter(4, 2*pi*EFilterBW,′s′); H_r =freqs(b,a,2*pi*Fgrid); ExpArg = −j*2*pi*Fgrid; ONE=ones(PtrnLength,1);%% Normalize the received OMA to 1. Estimate the OMA of the captured %%waveform by using a linear fit to estimate a pulse response, sythesize a%% square wave, and calculate the OMA of the sythesized square wave per%% Clause 52.9.5 ant=4; mem=40; %Anticipation and memory parameters forlinear fit X=zeros(ant+mem+1,PtrnLength); %Size data matrix for linearfit Y=zeros(OverSampleRate,PtrnLength); %Size observation matrix forlinear fit for ind=1:ant+mem+1 X(ind,:)=circshift(XmitData,ind-ant−1)′;%Wrap appropriately for lin fitend X=[X;ones(1,PtrnLength)]; %The all-ones row is included to computethe bias %Each column of Y is one bit period for ind=1:OverSampleRate Y(ind,:)=yout0([0:PtrnLength−1]*OverSampleRate+ind)′; endQmat=Y*X′*(X*X′){circumflex over ( )}(−1); % Qmat is the coefficientmatrix resulting from linear fit. Each column % (except the last) is onebit period of the pulse response. The last % column is the bias.SqWvPer=16; %Even number; period of the square wave used to compute theOMA SqWv=[zeros(SqWvPer/2,1);ones(SqWvPer/2,1)]; %One period of sqwave(col) X=zeros(ant+mem+1,SqWvPer); %Size data matrix for synthesisfor ind=1:ant+mem+1  X(ind,:)=circshift(SqWv,ind-ant−1)′; %Wrapappropriately for synthesis end X=[X;ones(1,SqWvPer)]; %Include thebaseline Y=Qmat*X;Y=Y(:); %Synthesize the modulated square wave, putinto one column %Set position in square wave to average overavgpos=[.4*SqWvPer/2*OverSampleRate:.6*SqWvPer/2*OverSampleRate];ZeroLevel=mean(Y(round(avgpos),:)); %Average over middle 20% of ″zero″run %Average over middle 20% of ″one″ run, compute OMAMeasuredOMA=mean(Y(round(SqWvPer/2*OverSampleRate+avgpos),:))−ZeroLevel;%% Subtract zero level and normalize OMA yout0 =(yout0-ZeroLevel)/MeasuredOMA; %% Compute the noise autocorrelationsequence at the output of the %% front-end antialiasing filter andrate-2/T sampler. Snn = N0/2 * fftshift(abs(H_r).{circumflex over( )}2) * 1/SymbolPeriod * OverSampleRate; Rnn = real(ifft(Snn)); Corr =Rnn(1:OverSampleRate/2:end); C = toeplitz(Corr(1:EqNf)); %% Compute theminimum slicer MSE and corresponding TWDP for the %% three stressorfibers % X matrix used in MSE calculation X = toeplitz(XmitData,[XmitData(1); XmitData(end:−1:2)]); Rxx = X′*X; %Used in MSE calculationChannelTWDP = [ ]; for ii=1:3 %index for stressor fiber  %% Propagatethe waveform through fiber ii.  %% The DC response of each fiber isnormalized to 1 (sum(PCoefs)=1)  PCoefs = FiberResp(ii+1,:)′;  Hsys =exp(ExpArg * Delays′) * PCoefs; Hx = fftshift(Hsys/sum(PCoefs));  yout =real(ifft(fft(yout0).*Hx));  %% Process signal through front-endantialiasing filter %%%%%%%%%%%%%%%  yout = real(ifft(fft(yout) .*fftshift(H_r)));  %% Compute MMSE-DFE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% The MMSE-DFE filter coefficients computed below minimize  %%mean-squared error at the slicer input. The derivation follows from  %%the fact that the slicer input over the period of the data sequence  %%can be expressed as Z = (R+N)*W − X*[0 B]′, where R and N are  %%Toeplitz matrices constructed from the signal and noise components,  %%respectively, at the sampled output of the antialiasing filter, W is  %%the feedforward filter, X is a Toeplitz matrix constructed from the  %%input data sequence, and B is the feedback filter. The computed W  %%and B minimize the mean squared error between the input to the  %%slicer and the transmitted sequence due to residual ISI and Gaussian  %%noise. Minimize MSE over 2/T sampling phase and equalizer delay and  %%determine BER  MseOpt = Inf;  % jj loop minimizes over sampling phase for jj= [0:OverSampleRate−1]−OverSampleRate/2   %% Sample at rate 2/Twith new phase (wrap around as required)   yout_2overT =yout(mod([1:OverSampleRate/2:TotLen]+jj−1,TotLen)+1);   Rout =toeplitz(yout_2overT, . . .    [yout_2overT(1);yout_2overT(end:−1:end-EqNf+2)]);   R = Rout(1:2:end, :);   RINV =inv([R′*R+PtrnLength*C R′*ONE;ONE′*R PtrnLength]);   R=[R ONE]; % Addall-ones column to compute optimal offset   Rxr = X′*R; Px_r = Rxx -Rxr*RINV*Rxr′;   %% kk loop minimizes MSE over equalizer delay,feedforward   %% coefficients, and feedback coefficients   for kk =1:length(EqDelVec)    EqDel = EqDelVec(kk);    SubRange =[EqDel+1:EqDel+EqNb+1];    SubRange = mod(SubRange−1, PtrnLength)+1;   P = Px_r(SubRange,SubRange);    P00 = P(1,1); P01 = P(1,2:end); P11 =P(2:end,2:end);    Mse = P00 - P01*inv(P11)*P01′;     if (Mse<MseOpt)    MseOpt = Mse;     B = −inv(P11)*P01′; % Feedback filter     XSel =X(:,SubRange);     W = RINV*R′*XSel*[1;B]; % Feedforward filter     Z =R*W - XSel*[0;B]; % Input to slicer     %% Compute BER usingsemi-analytic method %%%%%%%%%%%%%%%%%     MseGaussian =W(1:end−1)′*C*W(1:end−1);     Ber =mean(0.5*erfc((abs(Z-0.5)/sqrt(MseGaussian))/sqrt(2)));    end   end end  %% Compute equivalent SNR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %% Thisfunction computes the inverse of the Gaussian error probability  %%function. The built-in function erfcinv( ) is not sensitive enough  %%for the low probability of error cases.  if Ber>10{circumflex over( )}(−12) Q = sqrt(2)*erfinv(1-2*Ber);  elseif Ber>10{circumflex over( )}(−323) Q = 2.1143*(−1.0658−log10(Ber)).{circumflex over ( )}0.5024; else Q = inf;  end  %% Compute penalty%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  RefSNR = 10 * log10(Q0) + PAlloc; ChannelTWDP(ii) = RefSNR−10*log10(Q); end %% Pick highest value due tothe multiple fiber responses from ChannelTWDP. TWDP = max(ChannelTWDP)%% End of program

APPENDIX C SAMPLE MATLAB CODE, EXAMPLE 3 This MATLAB code illustrates anexample of various aspects of the invention. This code implements Method2 for automatically estimating OMA and baseline from the input waveformand includes temporal alignment of the synthesized square waveform anddata test pattern used for OMA calculation. function[xWDP,MeasuredxMA]=SFF8431xWDP(WaveformFile,EqNf,EqNb, . . .SymbolRate,Usage) %% Example calling syntax:%% [xWDP,MeasuredxMA]=SFF8431xWDP(′wavefile.txt′,14,5,10.3125,′Optical_WDP′) %% The fields in the example above should be replaced bythe actual values %% being used. WaveformFile should be thepath\filename for each waveform %% tested. The waveform consists of Nsamples per unit interval T, where %% N is the oversampling rate. Thewaveform is circularly shifted to align %% with the data sequence. Thefile format for the measured waveform is ASCII %% with a single columnof chronological numerical samples, in signal level, %% with no headersor footers. %% EqNf is the # of T/2-spaced feedforward equalizer taps;EqNb is the # of %% T-spaced feedback equalizer taps. %% SymbolRate isin gigabaud. %% Options for Usage are ′Optical_WDP′, ′Copper_WDP′, and′Copper_TWDP′. %% ′Optical WDP′ is used in support of clause 3 for%% measuring WDPi at the output of an optical TP3 tester, %% measuringWDPo at C′ of a linear optical module receiver, and %% calibrating WDPat C″ for testing a host for linear optical modules. %% ′Copper WDP′ isused in support of Annex E for %% measuring WDPi and calibrating WDP atB″ for a copper cable assembly, %% measuring WDPo at C′ of a coppercable assembly (C), and %% calibrating WDP at C″ for testing a host forcopper cable assemblies. %% ′Copper TWDP′ is used for measuring TWDP atB of a host for copper %% cable assemblies. %% Transmit data file: Thetransmit data sequence is 511 bit PRBS9 TWDP test %% patterns defined inTable 68-6. File format is ASCII with a single column %% ofchronological ones and zeros with no headers or footers. TxDataFile =′prbs9_950.txt′; %% Program constants %% OverSampleRate = 16; %Oversampling rate, must be even SymbolPeriod = 1/SymbolRate; % Symbolperiod is in ns Q0 = 7.03; % BER = 10{circumflex over ( )}(−12) %% Loadinput waveform and data sequence, generate filter and other matricesyout0 = load(WaveformFile); XmitData = load(TxDataFile); PtrnLength =length(XmitData); TotLen = PtrnLength*OverSampleRate; Fgrid =[−TotLen/2:TotLen/2−1].′/(PtrnLength*SymbolPeriod); %% Compute responseof 7.5 GHz 4th order Butterworth antialiasing filter a = [1 123.14077581.811 273453.7 4931335]; % Denominator polynomial b = 4931335; %Numerator for frequency response ExpArg = −j*2*pi*Fgrid; H_r =b./polyval(a,−ExpArg); %% Get usage parameters for the application[H_chan,Delays,PAlloc] = GetParams(Usage,ExpArg); N0 = SymbolPeriod/2 /(Q0 * 10{circumflex over ( )}(PAlloc/10)){circumflex over ( )}2; %% Setsearch range for equalizer delay, in symbol periods. Lower end %% ofrange is minimum channel delay. Upper end of range is the sum of the %%lengths of the FFE and channel. Round up and add 5 to account for the %%antialiasing filter. EqDelMin = floor(min(Delays)/SymbolPeriod);EqDelMax = ceil(EqNf/2 + max(Delays)/SymbolPeriod);ONE=ones(PtrnLength,1); %% Normalize the received xMA (OMA or VMA) to 1.Estimate xMA of the captured %% waveform by using a linear fit toestimate a pulse response, synthesize a %% square wave, and calculatethe xMA of the synthesized square wave per IEEE %% 802.3, clause 52.9.5.ant=4; mem=40; % Anticipation and memory parameters for linear fitX=zeros(ant+mem+1,PtrnLength); % Size data matrix for linear fitY=zeros(OverSampleRate,PtrnLength); % Size observation matrix for linearfit for ind=1:ant+mem+1  X(ind,:)=circshift(XmitData,ind-ant−1)′; % Wrapappropriately for lin fit end X=[X;ones(1,PtrnLength)]; % The all-onesrow is included to compute the bias for ind=1:OverSampleRate Y(ind,:)=yout0([0:PtrnLength−1]*OverSampleRate+ind)′; % Each column is1 bit end QmatY*X′*(X*X′){circumflex over ( )}(−1); % Coefficient matrixresulting from linear fit. Each %% column (except the last) is one bitperiod of the pulse response. The last %% column is the bias.SqWvPer=16; % Even number; sets the period of the sq wave used tocompute xMA SqWv=[zeros(SqWvPer/2,1);ones(SqWvPer/2,1)]; % One period ofsq wave (column) X=zeros(ant+mem+1,SqWvPer); % Size data matrix forsynthesis for ind=1:ant+mem+1  X(ind,:)=circshift(SqWv,ind-ant−1)′; %Wrap appropriately for synthesis end X=[X;ones(1,SqWvPer)]; % Includethe bias Y=Qmat*X;Y=Y(:); % Synthesize the modulated square wave, putinto one column Y=AlignY(Y,SqWvPer,OverSampleRate); %New subFunction inr1.2 avgpos=[0.4*SqWvPer/2*OverSampleRate:0.6*SqWvPer/2*OverSampleRate];ZeroLevel=mean(Y(round(avgpos),:)); % Average over middle 20% of ″zero″run % Average over middle 20% of ″one″ run, compute xMAMeasuredxMA=mean(Y(round(SqWvPer/2*OverSampleRate+avgpos),:))−ZeroLevel;%% Subtract zero level and normalize xMA yout0 =(yout0-ZeroLevel)/MeasuredxMA; %% Compute the noise autocorrelationsequence at the output of the front-end %% antialiasing filter andrate-2/T sampler. Snn = N0/2 * fftshift(abs(H_r).{circumflex over( )}2) * 1/SymbolPeriod * OverSampleRate; Rnn = real(ifft(Snn)); Corr =Rnn(1:OverSampleRate/2:end); C = toeplitz(Corr(1:EqNf)); %% Compute theminimum slicer MSE and corresponding xWDP X = toeplitz(XmitData,[XmitData(1); XmitData(end:−1:end+1−EqNb)]); Xtil =toeplitz(circshift(XmitData,EqDelMin), . . . XmitData(mod(−EqDelMin:−1:−(EqDelMax+EqNb),PtrnLength)+1)); Rxx = X′*X;% Used in MSE calculation %% Propagate the waveform through channel.yout = real(ifft(fft(yout0) .* fftshift(H chan))); %% Process signalthrough front-end antialiasing filter %%%%%%%%%%%%%%%%%% yout =real(ifft(fft(yout) .* fftshift(H_r))); %% Compute MMSE-DFE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% The MMSE-DFE filter coefficientscomputed below minimize mean-squared error %% at the slicer input. Thederivation follows from the fact that the slicer %% input over theperiod of data sequence can be expressed as Z = (R+N)*W − % X*[0 B]′,where R and N are Toeplitz matrices constructed from signal and %% noisecomponents, respectively, at the sampled output of the antialiasing %%filter, W is the feedforward filter, X is Toeplitz matrix constructedfrom %% the input data sequence, and B is the feedback filter. Thecomputed W and B %% minimize the mean square error between the input tothe slicer and the %% transmitted sequence due to residual ISI andGaussian noise. Minimize MSE %% over 2/T sampling phase and FFE delayand determine BER. MseOpt = Inf; for jj=[0:OverSampleRate−1]−OverSampleRate/2 % sampling phase  %% Sample atrate 2/T with new phase (wrap around as required)  yout_2overT =yout(mod([1:OverSampleRate/2:TotLen]+jj−1,TotLen)+1);  Rout =toeplitz(yout_2overT, [yout_2overT(1); yout_2overT(end:−1:end−EqNf+2)]);  R = Rout(1:2:end, :);  RINV = inv([R′*R+PtrnLength*CR′*ONE;ONE′*R PtrnLength]);  R=[R ONE]; % Add all-ones column to computeoptimal offset  Rxr = Xtil′*R; Px_r = Rxr*RINV*Rxr′;  %% Minimize MSEover equalizer delay  for kk = 1:EqDelMax−EqDelMin+1   SubRange =[kk:kk+EqNb];   SubRange = mod(SubRange−1, PtrnLength)+1;   P = Rxx -Px_r(SubRange,SubRange);   P00 = P(1,1); P01 = P(1,2:end); P11 =P(2:end,2:end);   Mse = P00 - P01*inv(P11)*P01′;   if (Mse<MseOpt)   MseOpt = Mse;    B = −inv(P11)*P01′; % Feedback filter    XSel =Xtil(:,SubRange);    W = RINV*R′*XSel*[1;B]; % Feedforward filter    Z =R*W - XSel*[0;B]; % Input to slicer    %% Compute BER usingsemi-analytic method %%%%%%%%%%%%%%%%%%    MseGaussian =W(1:end−1)′*C*W(1:end−1);    Ber =mean(0.5*erfc((abs(Z-0.5)/sqrt(MseGaussian))/sqrt(2)));   end  end end%% Compute equivalent SNR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Thisfunction computes the inverse of the Gaussian error probability %%function. The built-in function erfcinv( ) is not sensitive enough forlow %% probability of error cases. if Ber>10{circumflex over ( )}(−12) Q= sqrt(2) *erfinv(1-2*Ber); elseif Ber>10{circumflex over ( )}(−323) Q =2.1143*(−1.0658−log10(Ber)).{circumflex over ( )}0.5024; else Q = inf;end %% Compute penalty %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RefSNR = 10 *log10(Q0) + PAlloc; xWDP = RefSNR−10*log10(Q); %% End of main function%% GetParams subFunction function [H_chan,Delays,PAlloc] =GetParams(Usage,ExpArg);  switch Usage   case ′Optical_WDP′ % Identitychannel for optical    Delays = 0;    H_Chan = 1;    PAlloc = 6.5; %Total allocated dispersion penalty (dBo)   case ′Copper_WDP′ % Identitychannel for copper    Delays = 0;    H_Chan = 1;    PAlloc = 7.5;   case′Copper_TWDP′ % Cu TWDP stressor    ChanResp = [. . .     .0 .04849.09697 .14546 .19394 .24243 .29091 .33940 .38788, . . .     .43637.48485 .53334 .58182 .63031 .67879 .72728 .77576;     .0 .01500 .08200.18300 .22400 .16600 .09600 .06600 .04800, . . .     .03000 .02200.01900 .01700 .01200 .00800 .00800 .00700];    Delays = ChanResp(1,:);   PCoefs = ChanResp(2,:)′;    H_chan =exp(ExpArg*Delays)*PCoefs/sum(PCoefs); %With normalization    PAlloc =7.5;  end %% End of GetParams function %% AlignY subFunction function Y= AlignY(Y0,SqWvPer,OverSampleRate)  % Aligns the crossings of the xMAsquare waveform Y to the first sample  % in the vector, etc. The firstsample is equivalent to time = 0.  Y = Y0−mean(Y0); % AC-couple socrossings are at 0.  tmp = Y(1:end−1).*Y(2:end); % Crossings will bewhere this vector is <=0.  % Look only for the approximate crossing inthe middle by ignoring any  % within ~2 UI from the beginning. Due topossible misalignment of the  % captured waveform, this is the onlycrossing that is certain.  x =find(tmp(2*OverSampleRate:end)<=0,1)+2*OverSampleRate−1;  % Use linearinterpolation to find a more exact crossing point.  xinterp =x-Y(x)/(Y(x+1)−Y(x));  % Calculate the difference from the idealcrossing time then shift and  % re-sample to create the aligned squarewaveform.  xadj = xinterp−(SqWvPer/2*OverSampleRate+1);  Y =circshift(Y0,−floor(xadj));  Y = interpl(Y,(1:length(Y))′+xadj−floor(xadj),′spline′); %% End of AlignY function

What is claimed is:
 1. A method for testing a transmitter for an opticalfiber communications link, the method comprising: accessing a storedwaveform for the transmitter under test, wherein the stored waveformrepresents an optical output generated by the transmitter in response toa data test pattern applied to the transmitter; processing the waveformin software to simulate propagation of the waveform through a referencechannel and/or a reference receiver to generate a synthesized waveform;determining a temporal offset between logic transitions of thesynthesized waveform and the data test pattern; temporally shifting thesynthesized waveform to generate an aligned synthesized waveform basedon the temporal offset; computing an optical modulation amplitude (OMA)based on the aligned synthesized waveform and the data test pattern; andcalculating a performance metric for the transmitter based on the OMA.